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Simplification  of  Nonlinear  Indicial  Response  Models: 
Assessment  for  the  Two-Dimensional  Airfoil  Case 

Jerry  E.  Jenkins* 

Wright  Research  and  Development  Center,  Wright-Patterson  Air  Force  Base,  Ohio  45433 


Simplifications  to  the  functional  form  of  the  nonlinear  indicial  response  model  and  application  to  the  aero¬ 
dynamic  res^nse  due  to  arbitrary  motion  inputs  are  discussed.  Numerical  results  for  a  thin  two-dimensional  airfoil 
with  wake  distortion  nonlinearities  are  used  to  control  and  justify  the  assumptions  employed.  Useful  relationships 
between  indicial  response  parameters  and  steady-state  osciUatory  force  and  moment  data  are  obtained  by 
representing  the  indicial  response  with  a  Taylor  series  (in  terms  of  onset  motion  parameters)  and  approximating 
the  superposition  integral  with  an  asymptotic  expansion.  In  particular,  the  appearance  of  certain  harmonics  and 
their  variation  with  frequency,  amplitude,  and  mean  angle  of  attack  are  traced  to  specific  onset  parameters.  A 
nonlinear  parameter  identification  procedure  Is  proposed  by  which  active  onset  parameters  may  be  determined  from 
exi^rimental  data.  Restrictions  imposed  by  using  the  asymptotic  expansion  are  also  examined.  Nonlinear  stability 
derivatives  are  shown  to  be  related  to  specific  indicial  response  characteristics. 


Nomenclature 

A  =  amplitude  of  harmonic  oscillation  in  a,  rad 
=  section  lift  coefficient,  Whjqc 
=  apparent  mass  derivative,  1/rad 
=  mean  section  lift  coefficient  for  oscillatory  motion 
=  nonlinear  indicial  response,  section  lift  due  to  a 
c  =  chord  length,  ft 

F  =  deficiency  function,  difference  between  steady-state 
and  transient  indicial  responses,  I /rad 
Fq  =  linear  deficiency  function  evaluated  at  a  =  0,  1/rad 
/s  =  nonlinear  deficiency  function,  d^Fjda^  evaluated 
at  a  =  0,  1/rad^ 

,Vo  =  superposition  integral  linear  onset  parameter,  =  a 
=  superposition  integral  nonlinear  onset  parameter, 

=  a“a 

k  =  reduced  frequency,  ^wcIlU^ 

q  =  dynamic  pressure,  psf 

=  remainder  after  n  terms  of  partial  sum 
/  =  time,  nondimensionalized  by  2U^lc 

=  freestream  velocity,  ft/s 
a  =  angle  of  attack,  rad 

a.a  =  first  and  second  derivatives  of  a  with  respect  to  t, 

rad 

^0  =  mean  angle  of  attack  for  harmonic  motion,  rad 

t  =  nondimensional  time  at  step  onset 

Tj  =  elapsed  nondimensional  time  from  onset,  =  r  —  t 

o)  =  circular  frequency,  rad/s 

Introduction 

ARGE-AMPLITUDE  oscillatory  force  and  moment  data 
typically  exhibit  nonlinear  dependencies  on  amplitude 
and  frequency.  Often  these  effects  can  be  attributed  to  hys¬ 
teresis;  sometimes  “simple’’  nonlinear  variations  in  either  (or 
both)  static  or  dynamic  force  and  moment  data  are  the  cause. 
In  either  case,  aerodynamic  response  modeling  for  flight 
mechanics  analyses  is  a  difficult  problem. 
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The  current  work  focuses  on  establishing  relationships  be¬ 
tween  steady-state  oscillatory  force  and  moment  data  and 
models  capable  of  accounting  for  some  of  these  nonlinear 
effects.  Limits  of  applicability  for  these  models  is  also  ex¬ 
plored.  The  basis  for  this  study  is  nonlinear  aerodynamic 
modeling  work  reported  by  Tobak  and  Chapman  ‘  and  Tobak 
and  Schiff.- 

The  modeling  concepts  apply  equally  well  to  any  of  the  six 
force  and  moment  components  and  to  any  motion  input, 
including  control  deflection;  however,  lift  response  due  to 
angle  of  attack  is  used  throughout  to  illustrate  the  concepts. 

Nonlinear  Indicial  Responses 
The  nonlinear  indicial  response  is  the  most  generally  appli¬ 
cable  modeling  concept.  It  is  defined  in  terms  of  two  motions 
as  shown  in  the  following  sketch; 


oc 


where  a(^)  is  the  “reference  motion,”  defined  for  —  oo  <  /  ^  t; 
Ui  consists  of  a((;),  for  t  ^  t,  and  is  held  constant  at  a(T)  for 
t  >x;  and  consists  of  a((J),  for  /  ^  r,  but  jumps  instanta¬ 
neously  to  a(T)  -h  Aa  for  t  >  x.  The  nonlinear  indicial  lift 
response  is  the  limit,  as  step  height  Aa  approaches  zero,  of  the 
difference  between  the  corresponding  lift  time  histories. 

Evidence  supporting  the  need  for  such  a  definition  is  given 
in  Figs.  I  “3,  taken  from  Graham’s^  tow-tank  experiments 
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using  a  NACA  0015  airfoil.  The  point  is  that  motion  history 
clfects  can  have  a  profound  intlucncc  on  the  resulting  flow- 
fields,  even  when  there  arc  identical  onset  conditions.  In  fact, 
history  effects  are  completely  rcilecicd  in  the  wake  structure. 

Consider  two  angle-of-attack  time  histories,  as  shown  in 
Fig.  1,  which  are  interpreted  as  two  distinct  reference  mo¬ 
tions.  (Graham^  investigated  "poststall  maneuvers”  that  re¬ 
turned  to  zero  angle  of  attack  after  holding  at  high  a: 
however,  the  point  is  better  illustrated  by  examining  only  a 
portion  of  the  complete  motion.)  Flowficlds  at  the  termina¬ 
tion  of  each  maneuver  (x  =  90  deg)  arc  shown  in  Figs.  3.  The 
[low  is  left  to  nght  and  the  view  is  span  wise  from  the  root, 
slightly  above  and  behind  the  w’ing  section.  Thus,  only  the 
upper  surface  of  the  wing  section  is  visible;  the  leading  edge  is 
to  the  top  of  the  picture.  For  the  rapid  pitch-up  (alpha 
dot  =0.7),  the  “dynamic  stall”  vortex  is  less  diffuse,  more 
tightly  wrapped,  and  closer  to  the  airfoil  upper  surface  than 
for  alpha  dot  =  0.2.  Note  also  that  at  alpha  dot  =  0.7  the 
starting  vortex  associated  with  the  onset  of  the  pitch-up 
motion  (lower  right)  has  not  had  lime  to  be  swept  down¬ 
stream  and  can  still  exert  an  appreciable  inrluence  on  airfoil 
response.  Furthermore,  if  is  held  constant  at  this  point, 
these  distinctly  different  vortex  systems  will  also  exhibit 
unique  diffusion  and  convection  properties.  Correspondingly 
large  differences  in  lift  are  observed  (Fig.  2).  Clearly,  motion 
history  effects  can  be  important.  Therefore,  the  reference 
motion  is  required  to  establish  appropriate  flow  conditions 
from  which  to  measure  the  corresponding  indicial  response. 
The  ‘u  motion  then  establishes  ihe  duuv^e  in  response  due  lo  a 
small  periurbaiion  in  y.. 

As  shown  by  Tobak  and  Chapman,'  the  dependence  of 
nonlinear  indicial  responses  on  reference  motion  requires  that 
they  be  expressed  mathematically  as  a  functional,  i.e.. 


-  C^[yAt) 
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where  the  step  onset  is  at  t  =  r.  The  functional  dependency  on 
prior  motion  x(v;)  distinguishes  the  nonlinear  indicial  response 
from  its  linear  counterpart. 

Equation  ( 1 )  deiincs  the  Frechet  derivative  of  the  func¬ 
tional  Cr.[x,{9)],  as  noted  by  Tobak  and  Chapman.'  They  also 
suggest  that  bifurcations  of  physically  realizable  (asymptot¬ 
ically  stable  to  small  perturbations)  steady-state  solutions 
corresponding  to  x,  arc  signaled  by  loss  of  Frechet  differen¬ 
tiability.  Such  occurrences  are  of  considerable  interest  to  the 
study  of  hysteresis  effects. 

Possible  simplifications  to  the  functional  representation  for 
indicial  responses  have  been  suggested  by  Tobak  and  Schi.ff.' 
Suppose  that  the  motion  at;)  is  analytic  (in  the  strict  mathe¬ 
matical  sense)  over  the  interval  —  x  <  r  <  r.  In  this  case.  y{;) 
may  be  replaced  by  its  Taylor  series  expansion  about  =  r. 
Therefore. 

r),i(T),i{ :),...]  (2) 

where  the  independent  variables  :i(r),i(r),i'(r)....  are  the  co¬ 
efficients  of  the  Taylor  expansion.  On  physical  grounds,  the 
distant  past  is  expected  to  be  less  important  to  the  step 
response  than  motion  characteristics  just  prior  to  onset,  sug¬ 
gesting.  perhaps,  that  only  a  few'  Taylor  series  coetficients 
need  to  be  retained. 

Simpiihcations  based  on  unsteady  aerodynamic  characteris¬ 
tics  for  thin  tw'o-dimensiona!  airfoils  in  an  inviscid  and  incom- 
pres.sible  fluid  arc  studied  in  the  following  section.  This 
approach  allows  the  nonlinear  indicial  response  model  to  be 
examined  in  a  context  where  assumptions  may  be  controlled 
and  justified.  However,  results  reported  herein  are  not  appli¬ 
cable  to  the  extreme  conditions  evident  in  the  example  just 
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presented.  Notably  large-scale  separated  flows  with  corre¬ 
sponding  time-dependent  equilibrium  states  are  excluded. 

Nonlinear  thin-airfoil  characteristics  were  obtained  by  us¬ 
ing  NLWAKE,  a  computer  code  developed  by  Scott  and 
McCune.'^  NLWAKE  provides  numerical  solutions  for  a  non¬ 
linear  version  of  Wagner’s  integral  equation,  which  relates  the 
quasisteady  bound  vorticity  to  wake  vorticity.  To  do  this,  the 
shed  vorticity  is  discretized  and  the  wake  allowed  to  distort 
under  the  influence  of  both  bound  vorticity  and  other  wake 
elements.  Since  NLWAKE  disallows  separation,  only  nonlin¬ 
ear  wake  distortion  effects  are  included.  NLWAKE  cannot  be 
used  to  investigate  hysteresis  effects. 

Two-Dimensional  Airfoil  with  Wake  Distortion 

Nonlinear  step  responses  computed  with  NLWAKE  are 
shown  in  Figs.  4  and  5.  The  required  Frechet  derivatives  were 
computed  numerically  based  on  positive  and  negative  step 
heights  of  0.1  deg.  Figure  4  shows  indicial  responses  for 
constant  angle  of  attack  prior  to  onset.  In  Fig.  5,  a(<J)  is 
ramped  at  constant  rate  (positive  and  negative)  to  the  same 
onset  a  values.  Motion  history  effects  caused  by  wake  distor¬ 
tion  have  an  appreciable  influence  on  the  indicial  response. 
However,  history  effects  predicted  by  NLWAKE  correlate 
very  well  with  onset  angle  of  attack;  onset  rate  having  a 
negligible  influence  up  to  quite  large  values.  Thus,  for  the 
two-dimensional  airfoil  in  the  absence  of  separation,  indicial 
response  functionals  may  be  adequately  approximated  by  a 
function  only  of  a(T)  and  elapsed  time  from  onset  t  -  r.  This 
approximation  is  used  in  the  following  development  although 
separation  effects  are  probably  important  at  lower  onset  rates 
than  shown  in  Fig.  5. 

NLWAKE  cannot  predict  nonlinearities  in  the  static  lift- 
curve  slope.  However,  for  completeness,  this  possibility  is 
included,  although  coupling  between  quasisteady  characteris¬ 
tics  and  indicial  response  time  history  cannot  be  explored. 
Also,  there  are  situations  (e.g.,  following  a  Hopf  bifurcation) 
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Fig.  4  Onset  angle-of-attack  effect  on  step  response. 
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Fig.  5  Onset  rate  effect  on  indicial  response. 


where  the  response  is  time  dependent  even  in  the  limit  as 
elapsed  time  (since  step  initiation)  approaches  infinity.  How¬ 
ever,  a  steady-state  condition  is  assumed  to  exist.  In  these 
cases,  it  is  convenient  to  introduce  the  “deficiency  func- 
tion“  F\ 

/■[a(T),f  -t]  =  QJa(t),oo]  -  CiJpL{x),t  -  t] 

where  the  first  term  on  the  right  side  is  the  steady-state 
lift-curve  slope, 

QJaW.oo]  =  lim  QJa(T),/-T] 

r  —  x-*cc 

and  the  indicial  response  is  expressed  as  a  function  of  onset 
angle  of  attack  and  elapsed  time  as  suggested  above.  Thus  Eq. 
(2)  becomes 

QJa(c);/~T]-QJa(T),r-t] 

=  QJa(T),oo]  -  F[a(T),r  -  r]  (3) 

In  addition,  for  uncambered  airfoils  (consistent  with 
NLWAKE  restrictions),  both  the  deficiency  function  and 
static  lift-curve  slope  are  even  functions  of  a.  Furthermore, 
except  for  possible  bifurcation  points,  both  terms  on  the  right 
side  of  Eq.  (3)  are  expected  to  be  analytic.  Thus,  expanding 
them  in  Taylor  series  about  zero  of  angle  of  attack,  retaining 
only  even  powers  of  alpha,  Eq.  (3)  becomes 

ClMO'J  “  t]  -  QJa(T),r  -  r] 

=  Q,(0,oo)  -t-  0.5  [Q,(0,cx))]a=(T)  +  - 

00L‘ 

-Fo{0,t  -  T)  -  0.5F,(0,r  -  T)a=(T)  -h  -  (4) 

where 

f„(0,r-i)  =  F[a(r),r-T]U,.o 

and 

and  Fo(0,r— t)  are  the  classical  linear  terms  be¬ 
cause  they  are  independent  of  onset  conditions.  Note  that  if 
additional  onset  parameters  are  included,  multivariate  Taylor 
series  will  be  required  and  mixed  partial  derivatives  will 
appear. 

F2  was  computed  from  NLWAKE  step  responses  at  steady- 
state  onset  conditions  of  ±1.0  deg  and  is  shown  in  Fig.  6. 
Figure  4  shows  the  comparison  between  steps  initiated  at 
4.0-deg  alpha  as  calculated  directly  from  NLWAKE  and  as 
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Fig.  6  Alpha-squared  contribution  to  step  response. 
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given  by  Eq.  (4),  dropping  all  terms  higher  than  second  order 
and  retaining  only  the  linear  lift-curve  slope.  This  is  nearly  an 
exact  representation,  and  higher-order  terms  are  not  needed. 

When  Eq.  (4)  is  used  to  construct  lift  responses  to  arbitrary 
inputs,  as  discussed  in  the  following  section,  various  integrals 
of  Fq  and  will  be  required.  To  this  end,  it  is  useful  to  have 
analytical  models  for  each.  Exponential  approximations  for 
the  two-dimensional  airfoil  linear  deficiency  function  abound 
in  the  literature.  Jones’  expression,  as  given  by  Fung,^  is  used 
herein.  An  exponential  approximation  for  F,,  based  on  a 
least-squares  fit  to  NLWAKE  data,  is  shown  in  Fig.  6.  These 
approximations  are  of  the  form: 


(5a) 

F2^bi(\  (5b) 

where  =  1.037  ( 1/rad),  a.  =  0.0455,  Aj  =  0.300,  =  -148.3 

( 1/rad^),  b.  =  24.9,  and  b^  =  0.254. 

Equation  (4)  is  the  nonlinear  indicial  response  form  used 
throughout  the  remainder  of  the  paper;  principal  limitations 
on  its  use  are  the  following: 

1)  Reference  motions  must  be  analytic.  Indicial  response 
functionals  can  then  be  represented  as  functions  of  onset 
motion  parameters. 

2)  The  indicial  responses  are  assumed  to  approach  a 
steady-state  condition  as  elapsed  time  since  onset  becomes 
large. 

3)  In  quantitative  results,  terms  above  second  order  in  a 
have  been  neglected  based  on  comparisons  with  NLWAKE 
output  (see  Fig.  4).  This  has  been  justified  only  for  attached- 
flow  conditions  and  results  are  implicitly  restricted  to  a  rela¬ 
tively  small  angle-of-attack  range. 

4)  The  deficiency  function  and  the  static  lift-curve  slope 
have  been  expanded  in  Taylor  series  about  alpha  =  0  and  a 
symmetric  airfoil  is  considered;  only  even  terms  are  retained. 
This  leads  to  a  simpler  presentation  and  clearly  shows  that  the 
classical  linear  theory  is  contained  as  a  special  case.  Both 
terms  can  be  expanded  about  nonzero  values  of  alpha,  if 
desired. 


Response  to  Arbitrary  Motion 
Superposition  Integral 

As  shown  by  Tobak  and  Chapman,^  responses  to  arbitrary 
motion  inputs  can  be  calculated  by  using  a  generalized  super¬ 
position  integral.  If  a  bifurcation  in  the  steady-state  response 
occurs  at  /  ==  t^,  the  integral  has  the  form 


=  Q[M(0)]  -I-  lim 


da 

Q,Wc);/,t]  — dr 
lo 


da 


QJa(0;LT]  —  dr  +  AQ(;,a,) 
e  dr 


(6) 


w’here 


AQ(/;aJ  =  CJa(<);/,r,  +  e]  -  Q[a(c);/,r,  -  e] 

Thus,  the  superposition  integral  is  split  to  allow  the  solu¬ 
tion  to  change  discretely  to  a  new  equilibrium  state  (and  to 
avoid  the  singularity).  Since  NLWAKE  cannot  support  an 
investigation  of  hysteresis  effects,  it  is  assumed  that  no  bifur¬ 
cations  take  place  over  the  time  interval  0  to  t.  In  this  case, 
integration  can  proceed  directly  from  r  =  0  to  r  —  r  and 
AQ=0. 

Simplified  Forms 

Approximating  the  nonlinear  functional  by  Eq.  (4)  (drop¬ 
ping  terms  above  second  order)  and  changing  the  variable  of 
integration  to  elapsed  time  from  onset  ri  =  /  — r,  Eq.  (6) 


becomes 


Q(0  =  Q[r,a(0)I  +  Q,(0,oc)[a(r)  -  a(0)] 
a^(0  -a^ 


t\o)  r 

- ^o(0,t,)go(/  -r,)  dr, 

jo 


-0.5 


^:(0,ri)^2(^-ri)dr, 


(7) 


where 


(7a  “ 


If  Eq.  (7)  is  used  for  dynamic  analyses,  the  equations  of 
motion  become  integro-differential  equations  because  go  and 
g2  appear  under  the  integrals.  Further  simplification  would  be 
needed  to  avoid  this  complication.  For  linear  systems,  Etkin^ 
proposed  using  “aerodynamic  transfer  functions”  in  conjunc¬ 
tion  with  the  Laplace  transform  of  the  equations  of  motion. 
One  possible  approach  to  the  nonlinear  problem  is  to  repre¬ 
sent  Eq.  (7)  in  the  frequency  domain  by  using  higher-order 
Laplace  transforms  and  George’s  association  of  variables 
technique  (see  Ref.  7).  There  would  be,  then,  a  potentially 
manageable  algebraic  problem  with  advantages  similar  to 
Etkin’s  linear  system  approach. 

The  objective  here  is  limited  to  developing  relationships 
between  indicial  response  parameters  and  the  response  to 
arbitrary,  but  prescribed,  motion.  To  this  end,  the  integrals 

-Ti)  dT,;  1=0,2 

may  be  integrated  recursively  by  parts  to  give 
^.(T|)^,(r  -  T,)  dr,  =gi(0)laO) -giiOfniQ) 

+  imrA0  -  imAO)  +  -  -  [?,(0]/J0)  + 

(8a) 


where 


lAt)  = 

[■■■J 

r 

7\(Ti)dT,  -  dr, 

n  times 

/.(0)=J 

F;(T,)dT,  dt, 

n  times 

^|>  =  4(T|)^(^,(/  -T,)]  dr, 


(8b) 

(8c) 

(8d) 


Differentiability  of^j,  and  ^2  imposes  no  new  restrictions  since 
Eq.  (2)  requires  a(0  to  have  derivatives  of  all  orders. 

Finally,  a  series  representation  for  the  lift  response  to  an 
arbitrary,  but  specified,  motion  input  (in  the  absence  of 
bifurcation)  is  given  by  substituting  Eq.  (8a)  back  into  Eq. 
(7): 


i(t)  +  Q,(0,oo)a(0  + 


o^\t) 

6 


(9) 


where  a  steady-state  initial  condition  has  been  assumed  and 


I„  =  Ia„  and  J„  = 


sr 
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Note  that  an  apparent  mass  term  must  be  in¬ 

cluded  to  account  for  the  noncirculatory  part  of  the  lift 
response.  The  incompressible  apparent  mass  reaction  is  in¬ 
stantaneous;  therefore,  the  indicia!  response  requires  an  im¬ 
pulse  at  step  onset.  Neither  Fq  nor  is  includes  the  impulse. 
However,  it  can  be  included  as  a  Dirac  delta  function  embed¬ 
ded  in  the  indicial  response  representation.®  Subsequent  inte¬ 
gration  across  the  delta  function  produces  an  apparent  mass 
term  proportional  to  the  instantaneous  value  of  alpha  dot.  In 
the  linear  compressible  case,  Leishman^  gives  exponential 
models  for  the  noncirculatory  terms. 

Note  also  on  the  right  side  of  Eq,  (9)  the  top  line  (apparent 
mass  excluded)  gives  the  quasisteady  response.  The  last  two 
lines  account  for  the  transition  from  the  initial  quasisteady 
state  to  the  steady-state  condition;  i.e.,  I„(r)  and  J„(0  describe 
the  transient  response  and  vanish  for  large  t.  I„(0)  and  J„(0) 
provide  the  dynamic  response. 

If  a  bifurcation  in  possible  steady-state  solutions  occurs 
inside  the  time  interval,  the  superposition  integral  must  be 
split  in  accordance  with  Eq  (6).  The  integration  procedure 
outlined  earlier  remains  valid;  (he  impact  will  be  felt  in  two 
ways.  First,  by  definition,  static  lift-curve  slope  terms  will 
differ  on  each  side  of  the  bifurcation  (corresponding  changes 
in  the  deficiency  functions  are  expected).  Second,  because 
appears  in  the  limits  of  integration  (and  in  AQ),  terms  like 
giU  -  appear.  Since,  for  arbitrary  motion,  is 

not  known  a  priori,  the  time-domain  solution  becomes  awk¬ 
ward  for  flight  mechanists.  Aerodynamic  reaction  models,  as 
proposed  by  Hanflf,*®  avoid  this  difficulty. 

Properties  of  the  Series  Representation  . 

Equation  (9)  has  the  distinct  advantage  of  separating  mo¬ 
tion  variables  from  indicial  function  characteristics.  More 
specifically,  the  integrals  are  now  independent  of  motion 
input;  i.e.,  only  multiple  integrals  of  Fq  and  F2  are  required. 
This  property  leads  directly  to  the  desired  relationships  be¬ 
tween  specific  indicial  response  onset  parameters,  stability 
derivatives,  and  steady-state  response  to  oscillatory  motion 
inputs.  However,  simplification  has  been  obtained  at  the  cost 
of  generality.  Properties  of  both  series  (linear  effects  given  in 
terms  of  and  nonlinear  effects  given  in  terms  of  J„)  are 
examined  later. 

First,  a  special  word  of  caution  is  in  order.  It  may  be 
tempting  to  use  Eq.  (9),  retaining  higher-order  terms  to 
improve  accuracy,  directly  in  dynamic  analyses.  Even  in  the 
linear  case,  this  can  lead  to  disastrously  false  results.  The 
reason  is  that  expansion  of  the  integrals  alters  the  nature  of 
the  characteristic  equation,  changing  it  from  transcendental  to 
polynomial.  Each  additional  term  contains  higher-order  time 
derivatives  and  therefore  increases  the  polynomial’s  degree. 
Inevitably,  extraneous  roots  are  introduced.  Spurious  roots 
can  occur  in  awkward  places  (well  into  the  right-half  plane, 
for  example)  even  if  the  approximation  is  quite  good  within 
the  radius  of  convergence.  For  a  discussion  of  this  problem 
see  Ref.  1 1  and  the  references  cited  therein.  To  reiterate,  Eq. 
(9)  is  introduced  simply  as  an  analytical  tool  for  studying 
relationships  among  indicial  response  characteristics  and  re¬ 
sponses  to  prescribed  motions  such  as  those  encountered  in 
wind-tunnel  testing. 

Both  series  in  Eq.  (9)  are  asymptotic  expansions  [of  the  left 
side  of  Eq.  (8a)]  with  respect  to  sequences  defined  by  taking 
successive  time  derivatives  of  go  and  g^.  This  implies  that  they 
are  valid  for  “sufficiently  slow”  motions,  as  will  be  shown  for 
the  special  case  of  harmonic  oscillations  about  a  constant 
mean  angle  of  attack.  The  motion  also  has  an  arbitrary  phase 
angle  relative  to  a  reference  signal: 

a(/)  -  ao  +  co%{kt  +  (p)  ( lOa) 

goU)  =  -Ak  sin(/:/  +  <t>)  (10b) 


g^(t)  =  sin[3(/cr  +  <^)]  +  .^2  sin[2(kt  -h  (p)] 

$\n{kr  +  (p)}  (10c) 

where 

5i  =  [/4(4a5H-^2)]/4, 

Now  consider  the  sequences  of  functions 
d" 

/  =  o,2  (11) 

where  n  =  0,1, 2,3 . 

The  right  side  of  Eq.  (8a)  is  an  asymptotic  expansion  if  gi^ 
is  an  asymptotic  sequence  and  is  u(g,-«_i);  that  is,  if 

lim^i^=0  and  lim-A;^  =  0  for  all  n  (12) 

k-*0  g-^  ^~^^gi,n~\ 

General  expressions  for  go_„  and  g,,  can  be  obtained  by 
putting  Eqs.  (I Ob)  and  (10c)  into  Eq.  (il).  Similarly,  and  J„ 
are  evaluated  by  repeated  integration  of  Eqs.  (5a)  and  (5b)  in 
accordance  with  Eqs.  (8b)  and  (8c).  Proof  that  the  series 
satisfy  Eq.  (12)  is  then  straightforward.  For  example,  when  n 
is  odd,  g,  „  _  ,  is 

A-"( 3" -  '5,  sin3Q  +  2"  '  ‘5,  sin2n  +  5,  sinO) 
and  the  dominant  terms  of  vary  with  frequency  as 

k''-^'cos3Q  k"*'cos2Cl  k"  ^ '  cosQ 
k-+fi  ’  k^+f2  ’  k^+fj 

where 

Cl  =  kt  (p 

Note  that  /j,  /j,  and  are  independent  of  k.  Since  Q 
approaches  (p  in  the  limit,  the  nonlinear  terms  are  easily 
shown  to  satisfy  Eq.  (12).  A  similar  argument  holds  for  even 
values  of  «, 

The  most  important  property  of  asymptotic  expansions  is 
that  their  partial  sums  have  an  error  of  the  same  order  as  the 
first  term  omitted.  For  all  values  of  Eq.  (12)  guarantees 
that  the  error  can  be  made  arbitrarily  small  as  k  approaches 
zero.  However,  the  quality  of  the  approximation  depends  on 
the  convergence  properties  of  the  series  for  fixed  k.  If  partial 
sums  decrease  initially,  useful  approximations  can  be  obtained 
even  if  the  series  ultimately  diverge.  Thus,  the  question  of 
establishing  a  frequency  limit  for  practical  applications  is  best 
answered  by  examining  the  behavior  of  the  first  few  terms. 

To  this  end,  steady-state  lift  response  characteristics  for  an 
airfoil  oscillating  about  zero  mean  angle  of  attack  were  com¬ 
puted.  Motion  variables  are  given  by  Eqs.  (lOa-c)  with 
ao  =  0.  There  can  be  no  response  at  the  first  harmonic  fre¬ 
quency  2k  since  ^2  =  0  in  this  case. 

Approximate  lift  responses  were  computed  from  Eq.  (9) 
with  I„(0)  and  J„(0)  defined  by  Eqs.  (5a),  (5b),  and  (8c).  Both 
I«(0  and  J„(/)  were  set  to  zero  since  only  the  steady-state 
solution  was  desired.  Corresponding  exact  steady-state  solu¬ 
tions  were  computed  directly  from  Eq.  (7)  also  using  Eqs. 
(5a)  and  (5b)  for  Fq  and  is,  respectively.  Nonlinearity  in  the 
quasisteady  lift  curve  was  neglected  in  both  cases. 

Partial  sums  of  the  series  expansions,  normalized  by  corre¬ 
sponding  exact  solutions,  are  presented  in  Figs.  7-9.  In 
keeping  with  traditional  dynamic  testing  practice,  in-phase 
(cosine)  and  out-of-phase  (sine)  components  are  presented 
separately  because  they  can  be  associated  with  static  and 
dynamic  stability  derivatives.  (Splitting  the  series  is  permissi¬ 
ble,  as  the  sum  of  two  asymptotic  expansions  is  also  asymp¬ 
totic,  provided  they  are  defined  with  respect  to  the  same 
sequence.) 

(p 
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TERMS  IN  PARTIAL  SUM  -  N 
Fig.  7  Linear  in-phase  series  convergence. 
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TERMS  IN  PARTIAL  SUM  -  N 
Fig.  8  Linear  out-of-phase  series  convergence. 
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TERMS  IN  PARTIAL  SUM  -  N 
Fig,  9  Second  harmonic  in-phase  series  convergence. 

Figure  7  shows  partial  sums  representing  the  linear  in-phase 
component  as  a  function  of  frequency  and  the  number  of 
terms  retained.  Only  the  static  lift-curve  slope  is  included  at 
iV  =  0,  while  l2g“o.!  is  added  for  A  =  1,  and  so  on.  As  deter¬ 
mined  by  the  ratio  test,  series  for  the  linear  components 
converge  for  k  <02  =  0.0455.  Thus,  characteristics  just  below, 
at,  and  just  above  the  convergence  frequency  are  illustrated. 
Asymptotic  behavior  is  evident;  accuracy  is  improved  at  any 
N  as  frequency  decreases.  At  frequencies  above  the  conver¬ 
gence  frequency,  however,  the  partial  sums  diverge  immedi¬ 
ately.  If  only  the  quasisteady  term  is  retained,  common 
practice  for  rigid-body  dynamics,  errors  on  the  order  of  10 
and  1 1  %  will  be  incurred  at  reduced  frequencies  of  0.045  and 
0.050,  respectively.  Finally,  the  static  lift-curve  slope  is  within 
6%  and  convergence  is  rapid  (to  about  2%  with  one  addi¬ 
tional  term)  at  the  lowest  frequency  of  0.03. 

Convergence  properties  for  the  linear  out-of-phase  compo¬ 
nents  are  shown  in  Fig.  8  and  are  quite  similar  to  the  linear 
in-phase  component.  Greater  percentage  errors  for  the  initial 
term  (31,  64,  and  74%  in  order  of  increasing  frequency)  are 
expected  since  there  is  no  quasisteady  term  common  to  both 
exact  and  approximate  results. 


Properties  for  the  nonlinear  in-phase  contribution  at  3A:  are 
displayed  in  Fig.  9.  In  this  case,  divergence  occurs  at  frequen¬ 
cies  above  Ik  =  by  =  0.254.  Therefore,  computations  were  car¬ 
ried  out  at  reduced  frequencies  of  0.050.  0.0847,  and  0.090  to 
include  the  highest  frequency  shown  for  the  linear  case  and  to 
overlap  the  convergence  frequency.  Again,  divergence  is  im¬ 
mediate  outside  the  radius  of  convergence.  At  /c  ==  0.050  con¬ 
vergence  is  relatively  fast,  to  within  about  4%  in  three  terms; 
however,  the  initial  error  is  surprisingly  high  (35%)  consider¬ 
ing  that  this  frequency  is  well  removed  from  the  divergence 
boundary  of  0.0847. 

Similar  results  could  be  shown  for  the  remaining  nonlinear 
contributions;  but  the  point  is  that  the  series  diverge  for 
frequencies  greater  than  the  slowest  varying  exponential  in¬ 
volved  in  the  corresponding  indicial  response.  That  is,  if  for 
large  r,  the  indicial  response  varies  as  e'^°\  then  the  conver¬ 
gence  frequency  is  given  in  terms  of  the  indicial  response  time 
constant  by 


kT^  I 


where 


l/a 

For  responses  at  harmonics  of  the  forcing  function,  the 
harmonic  frequency  is  to  be  used. 


Relationship  to  Stability  Derivatives 

As  Elkin**  has  shown,  the  conventional  stability  derivative 
representation  can  be  derived  easily  from  the  premise  that  the 
aerodynamic  reactions  are  functionals  of  the  vehicle  state 
variables.  Two  assumptions  are  required: 

1)  The  motion  is  an  analytic  function  of  time. 

2)  Aerodynamic  reactions  are  analytic  functions  of  the 
instantaneous  state  variables  and  their  time  derivatives. 

The  first  assumption  allows  the  functional  to  be  replaced  by 
an  ordinary  function  of  the  stale  variables  and  their  deriva¬ 
tives;  the  second  permits  a  multivariate  Taylor  series  expan¬ 
sion  of  the  function  in  terms  of  its  arguments.  The  coefficients 
of  the  series  are  the  classic  stability  derivatives.  Only  the 
linear  terms  are  usually  retained  for  small  perturbation  analy¬ 
ses;  however,  nonlinear  derivatives  are  sometimes  included. 
Alternatively,  nonlinear  effects  can  be  included  by  treating  the 
derivatives  as  functions  of  the  variables. 

Similarities  between  Etkin’s  development  and  the  derivation 
of  Eq.  (4)  are  apparent.  The  significant  difference  is  that  the 
latter  is  an  approximation  to  a  functional  representing  the 
nonlinear  indicial  response  rather  than  one  representing  the 
total  response.  In  the  absence  of  bifurcation,  putting  Eq.  (4) 
into  the  superposition  integral  gives  the  total  response  to  an 
analytic  motion  input,  Eq.  (7).  In  this  section,  the  asymptotic 
approximation  to  Eq.  (7),  given  by  Eq.  (9).  is  used  to  briefly 
examine  the  implicit  limits  of  the  stability  derivative  ap¬ 
proach.  First,  to  gain  physical  insight,  Eq.  (9)  is  used  to 
replicate  a  previous  result  derived  by  Tobak  and  SchiflP  for 
the  case  of  slowly  varying  harmonic  oscillations  about  con¬ 
stant  mean  values. 

Recall  that  the  infinite  series  in  Eq.  (9)  are  expansions  for 
the  left  side  of  Eq.  (8a)  (with  /  =  0  and  /  =  2,  respectively) 
and  that  I„(/)  and  J„(0  tend  to  zero  as  /  goes  to  infinity.  Thus, 
taking  a  sufficiently  large  r,  say  and  retaining  only  the 
leading  terms  (/i  =  0)  of  the  series: 


F[x{t  -  r,),Ti]a(/  -Ti)  dti 

Jo 

^-a(O[/,(0)+0.5y,(0)«=(O]  (13) 
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Comparing  the  last  two  lines  of  Eq.  (13),  the  integral  is 
dominated  by  the  region  near  the  lower  limit  (t,  =  0).  This  is 
to  be  expected  on  physical  grounds  since  small  T;  corresponds 
to  the  recent  past  and  the  asymptotic  nature  of  the  series 
requires  the  motion  to  be  slow  in  some  sense.  This  suggests 
that  in  the  integrand  tj  can  be  neglected  compared  to  i.e., 


Consider  again  a  harmonic  oscillation  about  a  constant 
mean  angle  of  attack,  given  by  Eq.  (lOa).  Substituting  Eqs. 
(10b)  and  (10c)  into  Eq.  (9),  the  steady-state  lift  response  has 
the  form 

=  Gj  cosD  -f  G2  cos2fi  -h  G3  cos3Q 


-Ti),T,](i(f„  -T,)  di, 

Jo 


-  3i(c) 


^W^),T,]  dr, 


(14) 


Finally,  recalling  that  I,(/^)  and  Ji(r„)  are  negligible  and 

FMO.x,]  =  Fo(0,ti)  +0.5F2(0,T,)a^(O 

the  last  lines  of  Eqs.  (13)  and  (14)  are  easily  shown  to  be 
identical. 

Thus,  to  first  order,  lift  due  to  alpha  dot  is  proportional  to 
the  area  under  the  deficiency  function  evaluated  at  the  instan¬ 
taneous  angle  of  attack,  and  thus 


+  Hi  sinH  +  H2  sin2fi  -h  H3  sin3D  (16) 

The  time-invariant  response  arises  from  nonzero  values  of 
ao  and  is  simply  the  average  of  the  maximum  and  minimum 
quasisteady  lift  values; 

Qo  =  ^o{[(2a5  +  3A‘)/\2]C^^^^  +  Q,(0,oo)} 

From  Eq.  (9),  it  can  be  shown  that  all  of  the  in-phase 
contributions,  Gj  through  G3,  are  power  series  in  k  and 
contain  a  frequency  independent  quasisteady  term.  For  exam¬ 
ple, 

G,  =  ^  {Q,(0,oc.)  -  Uk^  +  Uk^  --  •••} 

+  +  (17a) 


ClM‘)]  =  dr,  =  1,(0)  +  0.5J,(0)a^(/J 

In  keeping  with  the  spirit  of  the  Taylor  series  interpretation  of 
the  stability  derivative  concept,  the  following  definitions  are 
perhaps  preferable: 


G2  =  -  AJ.k-  \6J^k*  (17b) 

where  and  are  evaluated  at  ti  =0. 

The  out-of-phase  components  are  similar  except  that  the 
leading  terms  are  first  order  in  frequency  and,  as  expected, 
contain  no  quasisteady  terms;  e.g., 


Q,  =  li(0) 


H2=  -(52/2)A:{Ji-4J3A’*-h-}  (I7c) 


Q,„=0.5:^^?^-0.5Ji(0) 

<7at7a- 


Of  course,  additional  information  is  obtained  by  including 
higher-order  terms  in  the  series  expansions.  For  example,  the 
second-order  terms  reveal  that  alpha-double-dot  effects  are 
related  to  the  second  integral  of  the  deficiency  function  again 
evaluated  at  a(r^),  as  shown  by  Eqs.  (15); 


Q,  =  l2(0)  and  Q,„=0.5J2(0)  (15) 

However,  both  an  additional  term  and  an  unexpected  identity 
appear; 


2Q,,.  =  Q,,,=0.5J,(0) 

Thus,  dynamic  stability  derivatives  and  the  dynamic  terms 
of  the  asymptotic  expansions,  terms  involving  I„(0)  and  J„(0), 
are  intimately  related.  When  the  asymptotic  approximation  to 
the  superposition  integral  breaks  down,  the  stability  deriva¬ 
tive  concept  does  also.  For  harmonic  motion  and  deficiency 
function  forms  given  by  Eqs.  (5a)  and  (5b),  this  occurs  at 
frequencies  at  or  above  the  largest  indicial  function  time 
constant  as  shown  in  the  preceding  section.  When  nonlineari¬ 
ties  introduce  harmonics,  the  harmonic  frequency  should  be 
compared  to  the  time  constant  for  the  corresponding  non¬ 
linear  component  of  the  deficiency  function. 


Relarionship  to  Steady-State  Oscillatory  Data 

Direct  measurement  of  nonlinear  indicial  responses  is 
difficult  at  best.  Furthermore,  suitable  parameter  identifica¬ 
tion  techniques  to  extract  them  from  experimental  data  have 
not  been  discussed  in  the  literature.  However,  an  assessment 
of  the  usefulness  of  approximating  the  indicial-response  func¬ 
tional  in  terms  of  onset  parameters  [Eq.  (4)]  needs  to  be 
accomplished  for  more  general  cases  than  studied  here.  The 
limited  objective  of  identifying  significant  onset  parameters 
from  steady-state  oscillatory  data  is  discussed  in  this  section. 


Now,  suppose  that  additional  terms  are  needed  in  Eq.  (4) 
to  represent  adequately  the  indicial  response.  For  example,  if 
onset  rate  is  important  and  its  effect  varies  linearly  with  angle 
of  attack,  terms  proportional  to  alpha  dot  and  the  product  of 
alpha  and  alpha  dot  should  be  considered.  That  is,  two 
additional  series  are  introduced  in  Eq.  (9)  involving  the 
motion  variables: 

and 

From  Eq.  (10a),  contributes  a  time-invariant  term  and  a 
first  harmonic.  Similarly,  contributes  a  constant  term  and 
components  at  the  fundamental,  first,  and  second  harmonic 
frequencies.  Thus,  including  these  terms  leaves  the  form  of 
Eq.  (16)  undisturbed;  however,  its  coefficients  (Gi,Hi,  etc.) 
then  contain  contributions  from  each  of  the  g,.  As  shown  in 
Eqs.  (17a-c),  the  effects  of  A  and  ao  are  separable  from 
frequency  effects;  i.e.,  they  may  written  in  the  form 

Gy  =  Go,yCo,y  +  G2.yC2,y  +  4"  ^4.y^4,y  (1^) 


where 


=  G,.  y(^,ao),  c,.  y  =  Cij{k) ;  j=\  ,2,3 

and  /  corresponds  to  onset  parameter  g/. 

Individual  effects  due  to  each  of  the  three  nonlinear  g/  on 
mean  lift  C/,y  and  /// y  are  summarized  in  Table  1.  Note  that 
each  onset  parameter  has  a  unique  signature  and  a  test  matrix 
can  be  designed  to  accentuate  the  differences. 


Table  1  Influence  of  onset  parameters 


Coefficient 

Sz 

8i 

^4 

0 

A^ 

A\ 

G,M, 

+  4A  ‘ag 

0 

A^ 

C,,//, 

A\ 

A^‘ 

A\ 

A^ 

0 

A^ 

138 


J.  E.  JENKINS 


J.  AIRCRAFT 


Given  a  candidate  set  of  onset  parameters,  the  identifica¬ 
tion  problem  consists  of  solving  for  the  unknown  c,  ^  on  the 
right  side  of  Eq.  (18).  Thus,  a  relatively  simple  harmonic 
analysis  of  data  obtained  by  conventional  techniques  is  suffi¬ 
cient  to  identify  active  onset  parameters.  Finally,  note  that 
this  analysis  requires  no  assumptions  about  the  form  of 
deficiency-function  time  history.  It  does  assume,  however,  that 
no  bifurcations  occur  within  the  test  range  and  requires  that 
the  series  in  Eq.  (9)  converge. 

Concluding  Remarks 

Earlier,  the  nonlinear  indicial  response  was  introduced  by 
Tobak,  Chapman,  and  Schiff  as  a  means  of  modeling  aero¬ 
dynamic  responses  in  nonlinear  flight  mechanics  problems. 
They  have  also  suggested  possible  simplifications  based  on 
relating  indicial  response  to  motion  variables  (and  their 
derivatives)  evaluated  at  step  onset. 

In  this  paper,  the  implications  of  using  a  Taylor  scries 
expansion  (with  respect  to  onset  parameters)  of  the  deficiency 
function  are  examined.  A  particularly  simple  expansion  is 
shown  to  be  applicable  to  thin  two-dimensional  airfoils  with 
wake  distortion  nonlinearities.  Furthermore,  given  the  series 
approximation  to  the  indicial  response,  the  generalized  super¬ 
position  integral  can  be  approximated  by  an  asymptotic  ex¬ 
pansion,  valid  for  sufficiently  slow  motions.  For  harmonic 
motion,  the  expansions  diverge  for  frequencies  greater  than 
the  slowest  varying  exponential  involved  in  the  corresponding 
deficiency  function.  The  combined  approximations,  for  the 
indicial  response  function  and  for  the  superposition  integral, 
lead  directly  to  relationships  between  conventional  stability 
derivatives  and  indicial  response  characteristics.  In  the  ab¬ 
sence  of  bifurcation,  these  relationships  can  be  used  to  iden¬ 
tify  dominant  onset  parameters  from  steady-state  oscillatory 
force  and  moment  data. 

Further  work  is  needed,  especially  concerning  bifurcation 
of  steady-state  responses  to  step  inputs  due  to  their  impor¬ 
tance  to  aerodynamic  hysteresis.  The  usefulness  of  represent¬ 
ing  the  indicial  response  as  a  function  of  onset  motion 
parameters  depends  on  obtaining  sufficient  accuracy  with  a 


very  limited  number  of  terms.  Assessments  based  on  experi¬ 
mental  data  at  realistic  reduced  frequencies  should  be  made. 
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Dynamic  wind-tunnel  test  results  of  a  65-deg  swept  delta  wing  are  reviewed.  These  tests  involved  body- 
axis  roiling  motions  at  moderate  (15-  to  35-deg)  angles  of  attack  in  both  the  Institute  for  Aerospace 
Research  2  X  3  m  low-speed  wind  tunnel  and  the  7  X  10  ft  Subsonic  Aerodynamic  Research  Laboratory 
facility  at  Wright- Patterson  Air  Force  Base.  They  included  static,  forced  oscillation,  and  free-to-roll 
experiments  with  flow  visualization.  Multiple  trim  points  (attractors)  for  body-axis  rolling  motions  and 
other  unusual  dynamic  behavior  were  observed.  These  data  are  examined  in  light  of  the  nonlinear  indicial 
response  theory.  The  analysis  confirms  the  existence  of  critical  states  with  respect  to  roll  angle.  When 
these  singularities  are  encountered  in  a  dynamic  situation,  large  and  persistent  transients  are  induced. 
Conventional  means  of  representing  the  nonlinear  forces  and  moments  in  the  aircraft  equations  of  motion, 
notably  the  locally  linear  model,  are  shown  to  be  inadequate  for  these  cases.  Finally,  the  impact  of  these 
findings  on  dynamic  testing  techniques  is  discussed. 


Nomenclature 

b  =  wingspan,  ft 

Cl,  C„  =  body-axis  rolling  moment  and  pitching  moment 
coefficients,  nondimensionalized  with  respect  to 
qSb  and  qSc,  respectively 

Qdyn  =  dynamic  rolling  moment  coefficient,  C/(f)  —  C;^u( 
Qsui  =  static  rolling  moment  coefficient,  i.e.,  equilibrium 
Cl  corresponding  to  the  instantaneous  value  of  roll 
angle 

Cift)  =  nonlinear  rolling  moment  indicial  response,  due  to 
step  input  in  roll  angle 
c  =  mean  aerodynamic  chord,  ft 

k  =  reduced  frequency,  <ob/2U„ 

q  =  freestream  dynamic  pressure,  Ib/fri 

Re  -  Reynolds  number 

S  =  wing  area,  fri 

t  =  time,  s 

=  freestream  velocity,  ft/s 

Xvb  =  vortex  breakdown  position,  fraction  of  root  chord 
aft  of  wing  vertex 

f  =  running  time  variable  denoting  motion  history, 

—oo  ^  <  T,  s 

cr  =  total  angle  of  attack,  body-axis  inclination  with 
respect  to  C/«,  deg 
r  =  time  at  step  onset,  s 
<f>  =  body-axis  roll  angle,  deg 

<f>Q  =  initial  roll  angle  for  free-to-roll  experiments,  also 
mean  roll  angle  for  harmonic  motion  experiments 
<o  -  circular  frequency,  rad/s 
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Introduction 

YNAMIC  coupling  between  aircraft  motion  and  aerody¬ 
namic  forces  and  moments  acting  on  the  aircraft  is  at  the 
heart  of  stability  and  control.  Maintaining  sufficient  fidelity  in 
aerodynamic  models  for  these  interactions  has  become  an  in¬ 
creasingly  difficult  problem  in  the  face  of  flight  envelope 
expansion. 

A  theoretical  method  for  studying  the  nonlinear  aspects  of 
the  flight  dynamics  problem  has  been  under  development  by 
Tobak^  and  his  colleagues  since  the  1960s.  Their  initial  ap¬ 
proach^  introduced  two  important  new  concepts:  1)  a  nonlinear 
indicial  response  and  2)  a  generalized  superposition  integral. 
As  with  linear  indicial  response  methods,  the  idea  is  to  rep¬ 
resent  aerodynamic  responses  (force  or  moment)  due  to  arbi¬ 
trary  motion  inputs  as  a  summation  of  responses  to  a  series  of 
step  motions.  The  nonlinear  indicial  response,  as  opposed  to 
its  linear  counterpart,  accounts  for  changes  induced  by  the  mo¬ 
tion  history  leading  up  to  step  onset. 

Subsequently,  results  from  the  growing  body  of  nonlinear 
dynamical  system  theory  were  used  to  greatly  strengthen  the 
model.^*'*  The  key  idea  of  these  extensions  has  been  to  accom¬ 
modate  the  existence  of  critical  states,  i.e.,  specific  values  of 
the  motion  variables  where  discrete  changes  in  static  aerody¬ 
namic  behavior  occur.  These  are  singular  points  that  require 
special  handling  in  the  superposition  integral.  Critical  states 
are  important  because  potentially  large  and  persistent  transient 
aerodynamic  effects  can  be  anticipated  when  an  aircraft  en¬ 
counters  them  in  a  dynamic  situation. 

Truong  and  Tobak^  have  also  demonstrated  that,  for  static 
aerodynamic  characteristics  that  are  time  invariant,  the  nonlin¬ 
ear  indicial  response,  together  with  the  generalized  superpo¬ 
sition  integral,  can  be  derived  directly  from  the  Navier- Stokes 
equations.  Thus,  the  theory  has  a  sound  mathematical  basis 
and  captures  the  flow  physics  as  well.  There  is  still  much  work 
to  be  done,  especially  for  cases  involving  time-dependent  equi¬ 
librium  states.  Nevertheless,  the  theory  is  rich  in  its  ability  to 
represent  a  wide  range  of  important  aerodynamic  nonlinearities 
that  can  be  encountered  in  maneuvering  flight. 

Independently,  Hanff®  proposed  the  reaction  hypersurface 
model.  As  opposed  to  the  time-domain  indicid  response 
model,  the  hypersurface  is  expressed  in  terms  of  a  set  of  in- 
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dependent  variables  consisting  of  the  instantaneous  values  of 
'the  motion  variables  and  their  time  derivatives.  It  was,  at  its 
inception,  designed  to  be  experimentally  based  and  primarily 
intended  for  simulations  of  aircraft  motion  in  nonlinear  settings 
where  the  classical  stability  derivative  approach  breaks  down.’ 
More  recent  work*  has  been  aimed  at  establishing  the  theoret¬ 
ical  connection  between  the  reaction  hypersurface  and  nonlin¬ 
ear  indicial  response  models. 

Experiments  designed  to  study  either  of  these  mathematical 
models  demand  a  dynamic  test  capability  that  can  efficiently 
collect  the  necessary  nonlinear  and  time-dependent  data.  A 
large-amplitude  high-rate  roll  oscillation  system,’  developed  at 
the  Canadian  Institute  for  Aerospace  Research  (lAR),  meets 
this  requirement. 

Hanff  and  Jenkins^°  used  this  rig  to  study  the  roll  dynamics 
of  both  a  65-deg  delta  wing  and  an  80-65-deg  double-delta 
wing  at  the  lAR,  Their  experiments  produced  some  extremely 
interesting  results  that  require  further  explanation. 

The  65-deg  delta  wing  configuration  was  found  to  have  mul¬ 
tiple  stable  trim  points  in  roll  (depending  on  roll-axis  incli¬ 
nation)  as  reported  by  Hanff  and  Ericsson."  Locations  for 
these  attractors  are  repeated  here  in  Table  1 .  They  argue  (based 
on  analysis  of  static  rolling-moment  data  at  cr  =  30  deg)  that 
asymmetric  vortex  breakdown,  induced  by  differing  effective 
sweep  angles  on  each  wing  panel,  is  the  root  cause.  However, 
the  dynamic  behavior  observed  in  free-to-roll  experiments  is 
harder  to  explain,  although  Hanff  and  Huang‘S  have  shown  that 
the  instantaneous  loads  are  largely  driven  by  the  dynamics  of 
leading-edge  vortex  breakdown. 

In  free-to-roll  tests,  the  model  is  given  an  initial  roll  dis¬ 
placement,  then  released  by  disengaging  a  remotely  actuated 
clutch.  The  model  is  then  free  to  roll  about  its  body  axis, 
restrained  only  by  a  small  amount  of  bearing  friction  in  this 
degree  of  freedom.  The  inertia  of  the  moving  system  was  ad¬ 
justed  to  ensure  that  the  free-to-roll  responses  were  in  the  same 
frequency  range  as  the  dynamic  force  measurements  (about  7.7 
Hz,  see  Table  2).  Since  frequency  is  inversely  proportional  to 
the  reference  length  at  fixed  flight  speed,  corresponding  full- 
scale  vehicle  frequencies  are  realistic. 

Two  free-to-roll  time  histories  for  the  65-deg  configuration 
at  cr  =  30  deg,  plotted  in  the  phase  plane,  are  shown  in  Fig. 
1.  Note  that  the  trajectory  for  the  — 58.3-deg  release  angle 
(solid  curve)  finds  the  stable  equilibrium  point  at  about  0-deg 
roll,  while  the  53.1-deg  release  trims  at  about  21  deg.  Both 
trajectories  pass  quite  close  to  attractors  (21  and  0  deg,  re¬ 
spectively)  with  very  low  rates,  but  do  not  trim  there.  This 
behavior  was  highly  repeatable.  Furthermore,  the  trajectories 
intersect  at  several  points.  Similar  intersections  of  phase-plane 
trajectories  (for  wing-rock  motions)  have  been  observed  only 


Table  1  RoU 
attractor  locations 


Angles,  deg 

Sting 

RoU 

20 

0 

25 

±1.5 

30 

0,  ±21 

35 

±11 

40 

0 

when  vortex  breakdown  occurs  over  the  wing.”  Clearly,  some 
phenomenon,  not  explicitly  accounted  for  in  the  two-dimen¬ 
sional  phase-plane  representation,  affects  the  motion.  Persis¬ 
tent  motion  history  effects,  perhaps  related  to  vortex  break¬ 
down  dynamics,  that  require  more  than  a  knowledge  of  the 
instantaneous  roll  angle  and  roll  rate  are  a  strong  possibility. 

Finally,  forced-oscillation  motions  about  a  zero  mean  roll 
angle  produce  distinctly  different  rolling-moment  responses 
than  those  measured  for  motions  with  mean  roll  angles  of  7 
deg  or  greater.'®  An  analysis  of  the  static  and  dynamic  force 
data*  suggested  that  this  behavior  can  be  explained  by  the  ex¬ 
istence  of  critical  states.  Indeed,  large  transient  effects  (with 
the  magnitudes  of  the  static  and  dynamic  effects  being  of  the 
same  order)  were  noted  in  the  response  following  encounters 
with  the  suspected  critical  states.  These  transients  persist  for 
at  least  a  quarter  cycle  at  /:  =  0.08. 

Follow-on  tests,  designed  to  investigate  the  cause  of  the  be¬ 
havior  noted  earlier,  were  conducted  in  the  Subsonic  Aero¬ 
dynamic  Research  Laboratory  (SARL)  wind  tunnel  at  Wright - 
Patterson  Air  Force  Base.  An  analysis  of  the  static,  dynamic, 
free-to-roll,  and  flow  visualization  data  taken  during  these  ex¬ 
periments  is  presented  later.  Relevant  aspects  of  critical-state 
theory  are  discussed  and  data  confirming  the  existence  of  roll- 
motion  critical  states  are  presented.  The  impact  on  dynamic 
testing  techniques  (including  data  collection  and  reduction) 
and  application  to  the  simulation  and  analysis  of  the  flight 
dynamics  problem  are  addressed  in  this  article’s  final  sections. 

SARL  Experiments 

lAR’s  high-amplitude,  high-rate  roll  apparatus  was  also  used 
in  the  SARL  tests.  A  comprehensive  experimental  program 
involving  over  800  runs  was  conducted  using  the  65-deg  delta¬ 
wing  configuration  (Fig.  2).  The  test  matrix  is  summarized  in 
Table  2. 

Since  SARL  is  an  open-return  atmospheric  tunnel,  the  LAR 
Mach  number,  Reynolds  number,  and  reduced  frequencies 
could  not  be  matched  simultaneously.  However,  test  conditions 
were  chosen  to  match  those  at  the  LAR  as  closely  as  possible. 
In  addition,  the  model  support  systems  and  tunnel  cross-sec¬ 
tion  geometry  are  quite  different  in  the  two  tunnels.  Thus,  the 


Table  2  SARL  test  conditions 


Test  type 

RoU  offset,  deg 

Amplitude,  deg 

Frequency,  Hz 

Total  angle  of 
attack,  deg 

Static  force 

-70^70 

NA 

NA 

15,  30,  35 

Dynamic  force 

0-42 

5-40 

1.1,  4.4,  7.7 

15,  30,  35 

Free-to-roll 

-65-65 

NA 

7.7 

30,  35 

Flow  visualization 

0-42 

5-40 

0,  1.1,  4.4.  7.7 

30.  35 

Surface  pressure 

0-42 

5-40 

0.  1.1,  2.2.  4.4,  7.7 

15.  30,  35 

iz 


270 


JENKINS,  MYATT,  AND  HANFF 


-  22.83" - 

Fig.  2  Model  geometry. 


very  good  facility-to-facility  repeatability  eliminated  the  pos¬ 
sibility  of  significant  support  and  wall  interference  effects. 

Dynamic  force  and  moment  measurements  were  taken  with 
the  model  forced  in  constant  amplitude  harmonic  motion.  Data 
were  taken  at  4.4  and  7.7  Hz  (k  =  0.08  and  0.14)  to  match 
lAR  conditions.  These  results  closely  tracked  static  measure¬ 
ments  except  when  the  motion  passed  through  small  roll  an¬ 
gles.*  Thus,  in  the  S  ARL  experiments,  dynamic  tests  were  also 
conducted  at  1.1  Hz  (it  =  0.02)  to  determine  how  the  extremely 
large  dynamic  effects,  incurred  at  small  roll  angles,  approach 
quasisteady  behavior  at  low  reduced  frequency. 

Laser-sheet  flow  visualization  data  and  unsteady  surface 
pressure  measurements  were  also  taken.  Thus,  an  extensive 
data  set  was  created  that  allows  a  coordinated  study  of  vortex 
dynamics  (including  breakdown)  and  the  resulting  unsteady 
aerodynamic  forces  and  moments, 

Free-to-roll  experiments  were  also  repeated  in  the  SARL 
tunnel.  These  data  provide  an  independent  check  on  the  ac¬ 
curacy  of  the  dynamic  force  measurements  (and  the  mathe¬ 
matical  model  used  to  represent  them)  since  the  measured 
forces,  together  with  the  known  model  and  test-rig  inertia,  can 
be  used  to  predict  the  free-to-roll  motion. 

Critical  States,  Theoretical  Basis 

A  development  of  the  nonlinear  indicial  response  and  critical 
states^'^’"^  is  beyond  the  scope  of  this  article.  However,  some 
key  results,  pertinent  to  the  present  discussion,  are  sum¬ 
marized: 

1)  The  nonlinear  indicial  response  (NIR)  is  represented 
mathematically  as  a  functional  to  incorporate  the  motion  his¬ 
tory  effect. 

2)  The  NIR  is  a  derivative,  called  the  Fr^chet  derivative,  of 
the  functional  representing  an  aerodynamic  response  in  terms 
of  its  motion  history.  It  is  the  limit,  as  step  height  goes  to  zero, 
of  the  incremental  response  (due  to  the  step  input)  divided  by 
step  height.  Following  Tobak*s  notation,*  the  rolling  moment 
due  to  an  infinitesimal  step  in  roll  angle  is  written 

QJit)  =  U  r] 

where  square  brackets  denote  a  functional;  the  first  argument 
is  the  independent  function  defining  the  motion  history  (roll 
angle  in  this  case);  and  arguments  following  the  semicolon 
give,  respectively,  the  times  at  which  1)  the  response  is  to  be 
evaluated  (observed)  and  2)  the  step  motion  was  initiated. 

Therefore,  the  function  <^(^)  is  to  be  interpreted  as  the  mo¬ 
tion  history  from  r  =  —  oo  to  step  onset  r,  and  the  motion  is  to 
be  held  constant  at  <^(t)  thereafter.  The  long-term  behavior  of 
the  aerodynamic  response  functional  (after  transients  due  to 
the  motion  have  died  away)  is  called  the  equilibrium  state. 

3)  If  the  Fr6chet  derivative  exists  everywhere  on  a  time  in¬ 
terval  (i.e.,  for  the  range  of  motion  variables  encountered  on 
that  interval),  the  generalized  superposition  integral  may  be 


used  to  construct  the  net  aerodynamic  response  (rolling  mo¬ 
ment  in  this  case)  over  the  interval.  Thus, 

C;(0  =  Q[<piO\  ^  0]  +  J  CJ(^(?);  r,  r]  —  dr 

Following  the  notation  introduced  previously,  the  first  term  on 
the  right-hand  side  (RHS)  is  the  rolling  moment  at  time  t  re¬ 
sulting  from  the  roll-angle  variation  which  is  the  motion 
history  prior  to  ^  =  r  =  0,  and  is  held  constant  at  for  all 
f  >  0.  The  functional  in  the  second  term  is  the  NIR,  as  defined 
earlier.  In  this  case,  t  is  the  variable  of  integration  and  the  time 
at  step  onset.  Thus,  the  integral  sums  the  effects  of  all  indicial 
responses  over  the  interval  0  to  t. 

4)  If,  on  the  other  hand,  there  are  specific  points  Tc  within 
the  interval  where  Fr^chet  differentiability  is  lost  (at  a  critical 
state  4>c),  the  integration  may  not  be  carried  beyond  the  instant 
at  which  a  critical  state  is  encountered  without  acknowledging 
the  existence  of  the  singularity. 

5)  Loss  of  Fr^chet  differentiability  is  handled  by  allowing 
the  equilibrium  response  to  change  discretely  to  a  new  state. 
Thus,  the  integral  must  be  split  to  isolate  the  critical  state,  i.e., 

r'"  d<^ 

Q(t)  =  C,[(^(f);  ^  0]  +  J  r,  r]  dr 

+  f  f.  t]  ^  dT  +  AC,(r,  4)  (1) 


where 

ACXr,  <l>c)  =  Q[<t>(i^);  r,  %  +  e]  -  r,  r,  -  e] 

(2) 

AC/,  as  given  by  Eq.  (2),  is  the  transient  response  associated 
with  4>c-  Note  that  it  depends  on  the  motion  history  from  — « 
to  just  beyond  Tc.  However,  its  effect  persists  for  times  t  >  r^. 

6)  Fr^chet  differentiability  may  be  lost  in  several  ways.*  A 
very  important  case  is  when  time-invariant  equilibrium  flows 
lose  their  analytic  dependence  on  a  motion  parameter  (body- 
axis  roll  angle  in  these  experiments).  There  are  at  least  two 
ways  this  can  happen,  both  involving  an  exchange  of  stability 
among  competing  equilibrium  flows: 

a)  Static  aerodynamic  responses  (i.e.,  responses  correspond¬ 
ing  to  solutions  of  the  time-invariant  form  of  the  Navier- 
Stokes  equations)  can  develop  a  fold  at  a  critical  value  of  the 
motion  parameter,  possibly  an  indirect  result  of  a  subcritical 
bifurcation.  The  response  slope  becomes  infinite  at  the  fold, 
invalidating  the  Frdchet  derivative.  Note  that  a  subcritical  bi¬ 
furcation  requires  the  existence  of  multiple  (nonunique) 
steady-state  solution  curves  on  both  sides  of  the  parameter’s 
critical  value. *^  However,  there  is  an  exchange  of  stability 
among  the  available  branches  at  the  critical  value.  The  initial 
or  basic  solution  becomes  unstable  when  the  parameter  is 
greater  than  the  critical  value.  Any  perturbation  will  force  a 
transition  to  a  new  (and  stable)  equilibrium  flow.  The  transition 
is  seen  as  a  discontinuous  jump  when  plotted  vs  the  motion 
parameter;  however,  its  time  histoiy  is  given  by  Eq.  (2). 

b)  Static  aerodynamic  responses  can  experience  a  supercrit¬ 
ical  bifurcation  at  a  critical  value  of  the  motion  parameter. 
Response-curve  slopes  become  discontinuous  at  such  points, 
again  invalidating  the  Fr^chet  derivative.  Nonunique  steady- 
state  solution  curves  are  required  in  this  case  also.*"*  The  ap¬ 
parent  discontinuous  jump  in  slope  is  caused  by  a  transition 
from  the  unstable  basic  solution  to  an  intersecting  (stable) 
equilibrium-solution  curve  as  the  motion  parameter  passes 
through  the  critical  value.  In  general,  the  intersecting  solution 
curves  are  smooth. 
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7)  Changes  in  flow  topology  (a  change  in  the  number  of 
singular  points  in  either  the  external  flow  or  in  skin-friction 
lines)  when  the  motion  parameter  reaches  a  critical  value  may 
signal  the  loss  of  Fr&het  differentiability.  Tobak  et  al.’  antic¬ 
ipate  that  equilibrium  responses  could  cease  to  be  analytically 
dependent  on  the  motion  parameter  at  such  points.  They  also 
state  that  there  could  also  be  “a  significant  increase  in  the  time 
required  for  the  , . .  response  to  reach  a  new  equilibrium 
state.” 

NIR  theory,  as  outlined  previously,  is  used  extensively  to 
interpret  the  new  experimental  results.  However,  there  are  im¬ 
portant  restrictions  implicit  in  the  analysis.  Static  force  data 
have  been  time  averaged  to  remove  the  effects  of  model  and/ 
or  sting  vibrations,  freestream  fluctuations,  etc.,  and  also  to 
stay  within  the  NIR  theoretical  framework  (which  in  its  current 
form  requires  time-independent  equilibrium  states).  Further¬ 
more,  the  relationship  between  time-averaged  force  and  mo¬ 
ment  behavior  and  topology  changes  in  the  corresponding 
mean  flow  deserves  closer  scrutiny  since  it  is  highly  desirable 
to  have  definitive  experimental  evidence  pointing  to  critical 
state  locations.  Therefore,  the  potential  impact  of  time-aver¬ 
aging  and  mean-flow  topology  changes  are  examined  in  the 
following  section.  This  is  done  by  using  a  well-known  flow  as 
an  example.  In  this  case,  the  parameter  controlling  the  flow- 
field  evolution  is  Reynolds  number  (as  opposed  to  an  evolution 
with  respect  to  a  motion  variable  as  discussed  earlier). 

Flow  About  a  Cylinder 

Consider  the  two-dimensional  flow  about  a  circular  cylinder 
at  Reynolds  numbers  ranging  from  zero  to  values  beyond  the 
onset  of  vortex  shedding.  As  Re  is  increased  from  zero,  the 
force  balance  (initially  between  only  pressure  and  viscous 
forces)  must  accommodate  the  increasing  effects  of  inertial 
forces.  At  very  low  Reynolds  numbers,  the  flow  remains  fully 
attached  to  the  cylinder  (as  shown  in  Fig.  3a);  the  flowfield 
topology  is  characterized  by  two  half-saddle  points’^  located 
at  the  fore  and  aft  stagnation  points  on  the  plane  of  symmetry. 
A  separation  bubble  appears  at  Re  ^  1  (based  on  the  cylinder 
diameter)  as  shown  in  Fig.  3b.  Associated  with  the  bubble  are 
three  half-saddle  points  on  the  surface  (an  increase  of  two) 
plus  a  saddle  point  in  the  downstream  flowfield  that  provides 
closure  for  the  bubble.  In  addition,  there  are  two  nodes  (de¬ 
generate  foci,  also  called  centers)  that  account  for  the  closed 
streamlines  within  the  bubble.  Thus,  both  flows  (with  and 
without  the  standing  eddies)  conform  to  the  topological  rule*^ 
for  the  flow  in  a  planar  slice  through  a  body: 


where  and  S5  denote  the  number  of  nodes  and  saddle 
points,  respectively,  in  the  flowfield  and  2/^  and  are  the 
number  of  half-nodes  and  half-saddle  points  on  the  body  sur¬ 
face,  respectively. 

Thus,  there  is  a  distinct  change  in  flowfield  topology  at  Re 
7  (going  from  two  topological  singularities  to  a  total  of 
seven,  including  three  off-surface  singularities).  However,  the 
changeover  is  perfectly  smooth.  There  are  no  competing  (non¬ 
unique)  equilibrium  solutions  for  1  <  Re  <  50;  i.e.,  the  topo¬ 
logical  change  does  not  involve  a  flowfield  instability  and  the 
possibility  of  a  supercritical  bifurcation  must  be  ruled  out 
However,  the  appearance  of  separated  flow  behind  the  cylinder 
implies  the  existence  of  reverse-flow  velocity  profiles  that  are 
more  susceptible  to  instability. 

At.  a  Re  ^  50,  the  flow  within  the  bubble  does  become 
unstable  and  vortex  shedding  begins.  The  resulting  equilibrium 
flow  is  time  periodic,  but  there  is  a  time-invariant  mean  flow. 
Numerous  ciculations  for  the  steady  flow  in  this  Re  range 
have  been  published. However,  experiments  by  Nishioka  and 
Sato*’  (Fig.  4)  clearly  show  that  the  separation  bubble  corre- 


Fig.  3  Steady  flow  about  a  cylinder:  a)  attached  flow,  Re  <1  and 
b)  with  separation  bubble,  Re  >  7. 
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Fig.  4  Separation  bubble  length  vs  Re  with  and  without  unsteady 
vortex  shedding. 


sponding  to  the  time-invariant  mean  flow  when  vortex  shed¬ 
ding  is  present  is  very  much  shorter  than  it  would  have  been 
in  the  absence  of  fluctuations.  Furthermore,  the  change  in  bub¬ 
ble  length  begins  abruptly  at  the  critical  Re  (corresponding  to 
the  onset  of  shedding),  and  entails  a  discontinuous  slope  with 
respect  to  Re.  A  corresponding  discontinuous  change  in  the 
slope  of  the  integrated  load,  i.e.,  the  mean  drag,  should  also 
be  expected. 

Note  that  there  has  been  no  change  in  the  number  and  types 
of  topological  singularities  present  in  the  mean  flow  at  Re 
50  (although  significant  changes  in  topology  are  evident  if  the 
flowfield  is  frozen  at  any  instant).  Thus,  two  time-invariant 
equilibrium  solutions  are  available  above  Re  50,  corre¬ 
sponding  to  the  cases  with  and  without  vortex  shedding. 
Therefore,  the  onset  of  a  time-periodic  equilibrium  condition 
is  seen  as  a  supercritical  bifurcation  of  the  mean  flow  wherein 
the  long  bubble  solution  loses  its  stability.  However,  there  is 
no  corresponding  change  in  mean-flow  topology. 

Recent  computations  by  Chen  et  al.*®  strongly  support  the 
notions  that  1)  there  are  no  bifurcations  associated  with  the 
appearance  of  the  separation  bubble  and  2)  loss  of  stability  in 
the  steady  flow,  through  a  Hopf  bifurcation,  is  associated  with 
the  onset  of  unsteady  vortex  shedding. 

Thus,  the  conclusions  are  as  follows: 

1)  Under  normal  conditions  (only  one  equilibrium  solution) 
a  change  in  mean-flow  topology  cannot  imply  a  loss  of  Frdchet 
differentiability;  i.e.,  there  is  no  corresponding  critical  state. 

2)  Supercritical  bifurcations  of  the  mean  flow  can  occur  (a 
critical  state)  without  a  coincident  change  in  mean-flow  to¬ 
pology  when  the  equilibrium  flow  becomes  time  dependent. 


Ill 
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3)  The  significance  of  a  topological  change  in  the  mean  flow 
is  that  if  inflexive  velocity  profiles  are  created,  then  the  re¬ 
sulting  mean  flow  is  more  susceptible  to  instability. 

These  observations  having  been  made,  the  discussion  now 
returns  to  the  rolling  delta  wing. 

Experimental  Evidence 

Based  on  the  theory  presented  previously,  aerodynamic  data 
in  the  presence  of  critical  states  should  exWbit  the  following 
traits: 

1)  Static  (fixed  vehicle  attitude)  flow  visualization  studies 
may  show  a  change  in  flow  structure  at  a  critical  state.  If  the 
visualization  technique  implicitly  involves  time  averaging 
(e.g.,  oil-flow  studies  of  skin-friction  lines),  a  topology  change 
is  neither  necessary  nor  sufficient  for  the  existence  of  a  critical 
state.  However,  a  mean-flow  topology  change  can  indicate  a 
susceptibility  to  instability  and  may  be  closely  associated  with 
a  subsequent  critical  state,  especially  if  flowfield  variables 
(such  as  pressure  gradient)  are  highly  sensitive  to  small 
changes  in  the  parameter. 

2)  Static  force  and  moment  data  should  exhibit  nonanalytic 
behavior  across  critical  states;  i.e.,  there  should  be  a  discon¬ 
tinuity  in  the  variation  of  the  force  data  and/or  their  derivatives 
with  respect  to  the  motion  variable  if  the  equilibrium  flow  is 
time  invariant.  The  discontinuity  is  located  at  the  critical  state. 
If  the  equilibrium  response  is  not  time  invariant,  a  disconti¬ 
nuity  may  still  be  apparent  in  the  mean  (time-averaged)  loads, 
as  in  the  example  presented  earlier.  This  latter  possibility  is  a 
likely  consequence  of  the  basic  flow  becoming  unstable, 
thereby  sending  the  unsteady  equilibrium  solution  on  a  new 
and  stable  path. 

3)  Transient  effects  should  be  observed  following  dynamic 
critical-state  encounters.  The  transient,  AC,  in  Eqs.  (1)  and  (2), 
will  in  general  depend  on  motion  history. 

SARL  data  (for  the  rolling  65-deg  delta  wing)  pertinent  to 
each  attribute  are  discussed  later.  All  data  discussed  here  were 
taken  at  the  same  condition  (0.3  Mach  number  and  cr  =  30 
deg).  For  these  discussions,  the  term  static  refers  to  the  model 
being  held  fixed  (within  normal  experimental  constraints)  with 
respect  to  the  freestream  and  does  not  imply  a  time-invariant 
equilibrium  state.  All  measurements  taken  under  static  condi¬ 
tions  are  time-averaged  unless  explicitly  stated  to  be  otherwise. 

Static  Flow  Structure 

Equilibrium  vortex  breakdown  locations  for  the  left  wing  as 
a  function  of  static  roll  angle,  from  Hanff  and  Huang, are 
shown  in  Fig.  5.  (For  positive  roll  angles,  the  left  wing  is  on 
the  lee  side;  i.e.,  it  has  been  rolled  away  from  the  freestream 
velocity  vector.  Conversely,  it  is  the  windward  wing  for  neg- 


Fig.  5  Time-averaged  vortex  breakdown  position:  left  wing 
panel,  <r  =  30  deg. 
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Fig.  6  Rolling-moment  coefficient  power  spectral  density:  a) 

=  20  and  b)  5  deg;  and  rolling  moment  coefficient  autocorrelation: 
c)  =  20  and  d)  5  deg. 


ative  roll  angles.)  The  solid  line  is  a  linear  least-square-error 
fit  to  the  circular  data  points  (correlation  coefficient  0.996). 
Note  that  the  experimental  point  at  ^  =  5  deg  (triangle)  departs 
significantly  from  the  linear  fit  Thus,  in  agreement  with  Hanff 
and  Huang’s  interpretation  of  Wentz’s'^  data,  breakdown  lo¬ 
cation  is  seen  to  be  a  nonlinear  (perhaps  discontinuous)  func¬ 
tion  of  static  roll  angle  in  the  4-  to  5-deg  range.  For  ^  greater 
than  5  deg  the  breakdown  point  is  well  aft  of  the  trailing  edge. 
In  addition,  vortex  breakdown  reaches  the  vertex  at  about  <f>  = 
—  13  deg,  as  suggested  by  the  linear  extrapolation.  When  the 
leading-edge  vortex  structure  on  both  wings  is  considered,  the 
conditions  |<^|  «  5  and  13  deg  are  candidate  critical  states. 

The  first,  \<i>\  5  deg,  may  be  a  critical  state  principally 

because  a  discontinuous  rearward  shift  of  the  breakdown  point 
on  the  lee  wing  would  cause  a  jump  (increase)  in  the  equilib¬ 
rium  lifting  pressures  on  that  side  (aft  of  about  80%  chord 
from  Fig.  5).  Any  discontinuities  in  the  forces  (or  their  deriv¬ 
atives)  with  respect  to  would,  of  course,  invalidate  the 
Fr6chet  derivative.  Time-averaged  static  force  data  supporting 
this  view  are  discussed  later. 

Secondly,  13  deg  is  almost  certainly  a  critical  state 

because,  if  the  breakdown  point  reaches  the  vertex,  the  stag¬ 
nation  point  in  the  vortex-core  axial  fiow^®  is  lost.  The  precise 
roll  angle  where  this  occurs  is  unknown  due  to  flow  visuali¬ 
zation  difficulties  that  require  an  extrapolation  of  the  data  to 
the  apex.  However,  the  breakdown  point  will  unquestionably 
reach  the  vertex  for  a  sufficiently  large  roll  angle. 

Static  Force  and  Moment  Behavior 

Unsteady  aerodynamic  force  and  moment  measurements, 
taken  under  static  conditions  where  vortex  breakdown  is  ab- 
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Fig.  7  a)  Static  (time-averaged)  roUing-moment  coefficient  data,  o-  =  30  deg;  b)  fuU-range  polynomial  fit  to  rolling-moment  data:  a  = 
30  deg,  |<^|  <  20  deg;  c)  time-averaged  pitching-moment  coefficient  data:  o-  =  30  deg,  |<^|  <  20  deg;  and  d)  polynomial  fits  for  the  ranges 
\<f>\  <  4.5  deg  and  |^|  >  4.5  deg. 


sent,  are  characterized  by  relatively  small,  broad-band,  random 
fluctuations  about  a  well-defined  mean  value.  (Mean  values 
calculated  for  overlapping  13-s  samples  repeat  quite  well.) 
Power-spectral  density  (PSD)  and  autocorrelation  plots  for  the 
rolling-moment  coefficient  for  such  a  case  (<7  =  30  deg,  4>  = 
20  deg)  are  presented  in  Figs.  6a  and  6c.  When  vortex  break¬ 
down  is  near  the  trailing  edge  of  the  left  wing  and  at  about 
40%  root  chord  on  the  right  wing,  o*  =  30  deg  and  =  5  deg, 
there  is  a  marked  increase  in  the  low-frequency  content  of  the 
PSD  (Fig.  6b),  i.e.,  for  nondimensional  frequencies,  (/  b)/ 
(2U^X  less  than  about  0.04.  As  a  result,  the  nns  value  of  the 
fluctuations  goes  from  0.0025  to  0.0033  (a  32%  increase)  and 
the  autocorrelation.  Fig.  6d,  takes  on  a  distinctly  quasiperiodic 
appearance. 

Time-averaged  static  rolling-moment  coefficient  vs  roll  an¬ 
gle  is  presented  in  Fig.  7a.  Also  shown  are  the  critical  states 
locations  corresponding  to  the  loss  of  a  coherent  vortex  on  the 
windward  wing  (breakdown  reaches  the  vertex).  The  degree 
of  unsteadiness,  in  terms  of  the  rms  value  of  the  fluctuations 
about  the  mean,  is  shown  by  error  bars  in  Fig.  7d. 

If  critical  states  exist  in  the  roll-angle  range  tested,  then  the 
function  representing  the  rolling  moment  (or  any  other  force/ 
moment)  cannot  be  analytic.  Given  only  values  for  a  function 
at  a  discrete  number  of  points,  it  is  impossible  to  prove  that  it 
is  (or  is  not)  an  analytic  function.  However,  a  cursory  inves¬ 
tigation  of  its  behavior  is  worthwhile. 

If  a  function  is  analytic  it  can  be  represented  by  a  Taylor 
series.  Therefore,  a  stepwise  regression  analysis  was  used  to 
look  for  a  polynomial  (least-squares  fit)  representing  C,  as  a 
function  of  since  the  polynomial  can  be  interpreted  as  a 
truncated  Taylor  series.  Thirty-two  terms  (each  an  odd  Legen¬ 
dre  polynomial^*)  were  included  in  the  list  of  possible  contrib¬ 


utors  to  the  regression  equation.  Test  data  points  were  inter¬ 
polated  (with  a  cubic  spline)  at  255  points  distributed  in  a 
geometric  sequence  beginning  at  the  origin  and  increasing  in 
either  direction;  i.e.,  the  ratio  of  to  <t>i  was  held  constant. 
The  proportionality  constant  between  successive  interpolation 
points  was  systematically  varied  such  that  intervals  between 
successive  points  (^,+,  —  0.)  at  the  extremes  ranged  from  2 
to  40  times  the  intervjd  at  the  origin  <f>i  -  This  was  done 
to  give  additional  weight  to  points  near  the  origin,  where  the 
rolling  moment  is  rapidly  varying,  while  keeping  enough  den¬ 
sity  at  the  extremes  to  prevent  the  polynomial  from  oscillating. 
Values  between  30-40  were  used  in  the  final  results  discussed 
later. 

The  solid  curve  (labeled  Full  Range  Fit)  in  Fig.  7a  is  the 
polynomial  fit  over  the  complete  roll-angle  range.  Thirteen 
terms  were  included  in  the  final  regression  equation;  the  rest 
were  rejected  because  they  offered  no  improvement  in  the  cor¬ 
relation.  Note  that  the  stable  trim  point  at  =  0  deg  was  not 
captured  despite  the  heavy  emphasis  given  to  points  in  this 
region  (160  of  255  points  were  within  |<^|  <  16  deg).  However, 
this  result  may  be  influenced  by  a  gap  in  test  data  coverage 
on  either  side  of  the  critical  state  at  about  13  deg  (raising 
questions  about  the  validity  of  the  spline  fit  in  this  region). 
Further  static  testing  is  required  to  determine  the  actual 
behavior. 

An  expanded  view  of  the  region  “20  deg  ^  <(>  ^  20  deg 
is  presented  in  Fig.  7b.  Also  shown  is  the  cubic  spline  fit  and 
the  possible  critical  state  at  5  deg  (from  flow  visualization 
results).  Note  the  separation  of  data  points  for  |<^|  <  4  deg 
and  for  5  deg  ^  |<5>|  ^7  deg.  Moreover,  Fig.  7c  (pitching 
moment  vs  roll  angle)  reinforces  the  notion  that  there  is  a 
discontinuity  between  the  two  groups.  Both  effects,  a  nose 
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Fig.  8  Rolling-moment  responses:  harmonic  motion,  a)  vertical  bar  denotes  critical  state,  b)  no  critical  state  encounters,  and  c)  vertical 
bars  denote  critical  states;  and  d)  dynamic  rolling-moment  responses:  harmonic  motion,  vertical  bars  denote  critical  states. 


down  pitching  moment  increment  and  a  positive  (right  wing 
down)  rolling-moment  increment  as  the  roll  angle  increases 
from  +4  to  +5  deg,  are  consistent  with  the  aft  movement  of 
the  vortex  breakdown  point  on  the  lee  wing  noted  earlier. 

Finally,  as  shown  in  Fig.  7d,  considerable  improvement  in 
the  polynomial  fit  is  possible  if  the  independent  variable  range 
is  broken  into  regions  and  distinct  polynomials  are  used  for 
each  of  them.  Interestingly,  if  breakpoints  are  included  only  at 
<l>  =  ±4.5  deg  to  model  apparent  discontinuities  there,  the 
polynomial’s  slope  at  4.5'^  deg  contrasts  sharply  with  test  data 
at  5,  6,  and  7  deg.  An  additional  breakpoint  for  the  critical 
state  at  approximately  13  deg  would  resolve  this  mismatch  in 
slopes.  In  summary,  the  inability  to  perform  polynomial  fits 
(without  breaking  the  interval  into  sm^er  pieces)  supports  the 
notion  that  critical  states  exist  near  the  locations  cited  previ¬ 
ously.  Again,  additional  static  data  are  needed  to  fully  justify 
this  conclusion. 

Thus,  for  |<^1  decreasing  through  5  deg  (breakdown  moving 
from  the  wake  to  a  position  over  the  wing),  static  flow  visu¬ 
alization  and  force/moment  data  suggest  that  a  change  in  flow 
topology  is  quickly  followed  by  a  subcritical  bifurcation  in¬ 
volving  an  oscillatory  instability.  (Chadeijian  and  Schiff’s^ 
computational  results  reinforce  the  view  that  there  are  large- 
scale  fluctuations  present  in  the  equilibrium  state.)  However, 
as  the  vortex  breakdown  approaches  the  wing  apex,  13 
deg,  the  bifurcation  that  occurs  would  appear  to  be  supercrit¬ 
ical. 

Dynamic  Forces  and  Moments 

If  the  conditions  outlined,  \<j>\  5  and  13  deg,  are  indeed 

critical  states,  there  may  be  discernible  transients  associated 
with  them.  A  previous  analysis*  of  lAR  data  suggested  that  at 
least  one  of  these  conditions  contributed  to  significant  dynamic 
effects.  In  this  section,  SARL  rolling-moment  measurements 
under  dynamic  conditions  are  examined  for  this  behavior. 


Recall  that  the  dynamic  data  are  steady-state  responses  to 
harmonic  motion,  i.e.,  starting  transients  have  dissipated  and 
the  aerodynamic  response  for  each  cycle  is  identical.  Thus,  if 
time  scales  appropriate  for  the  aerodynamic  response  are  suf¬ 
ficiently  large  compared  to  the  period  of  the  motion,  the  mea¬ 
sured  responses  represent  an  aggregate  of  effects  that  were 
initiated  during  earlier  cycles.  A  significant  frequency  effect 
will  be  observed  under  these  conditions.  Furthermore,  transient 
responses  to  single  events  occurring  at  discrete  points  during 
each  cycle  become  apparent  as  the  frequency  of  the  motion  is 
decreased. 

Dynamic  rolling-moment  data  taken  for  a  12-deg  amplitude 
body-axis  rolling  motion,  centered  about  a  mean  roll  angle  of 
14  deg,  is  presented  in  Fig.  8a.  The  abscissa  is  the  argument 
{(ot)  of  the  cosine  function  that  defines  the  motion.  Thus,  pre¬ 
cisely  one  cycle  of  motion  is  presented  regardless  of  frequency. 
Dynamic  data  taken  at  three  reduced  frequencies  (k  -  0.02, 
0.08,  and  0.14,  respectively),  are  shown.  In  addition,  the  roll¬ 
ing-motion  time  history  and  static  data  (plotted  as  functions  of 
the  instantaneous  roll  angle)  are  presented  for  reference. 

Note  the  dramatic  differences  in  waveform  between  re¬ 
sponses  at  the  two  highest  reduced  frequencies  and  the  k  = 
0.02  data  (where  distinct  transients  originating  at  critical  states 
become  more  apparent).  The  critical-state  encounter  at  an  (ot 
of  about  95  deg  (highlighted  by  the  vertical  bar)  is  readily 
discemable.  However,  later  in  the  cycle,  where  transient  effects 
overlap,  a  positive  identification  of  distinct  events  is  more 
difficult. 

Significantly,  responses  for  all  three  frequencies  follow  the 
static  data  closely  for  ot  in  the  range  0  to  about  95  deg  (26  deg 
>  2:  13  deg),  then  depart  from  the  static  curve.  In  addition, 

dynamic  responses  at  all  frequencies  approach  the  static  value 
at  the  end  of  the  cycle,  the  deviation  increasing  with  frequency. 
Thus,  the  rolling-moment  response  appears  to  be  essentially 
quasisteady  when  transients  due  to  critical-state  encounters  have 
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had  time  to  die  out.  This  notion  is  supported  by  the  data  pre¬ 
sented  in  Fig.  8b,  which  show  static  and  dynamic  rolling-mo¬ 
ment  data  for  12-deg  oscillations  about  a  mean  roll  angle  0o  of 
28  deg.  Roll  angles  traversed  by  this  motion,  16-40  deg.  pre¬ 
clude  any  encounters  with  suspected  critical  states.  Differences 
between  the  static  and  dynamic  data  are  likely  within  interpo¬ 
lation  errors  caused  by  the  sparsely  spaced  static  data  in  this 
region.  Again,  if  there  are  no  recent  critical-state  encounters,  the 
rolling-moment  behavior  is  quasisteady 

The  situation  shown  in  Fig.  8c  is  somewhat  different;  in  this 
case,  the  motion  (12-deg  amplitude  about  a  3-deg  offset)  is 
centered  between  critical  states.  Therefore,  the  order,  number, 
and  elapsed  times  between  encounters  are  different  than  those 
in  Fig.  8a.  Also,  roll  rates  at  |<^|  =  5  deg  are  much  higher  than 
before  (the  converse  is  true  at  ^  =  12  deg).  Definite  breaks  in 
the  rolling-moment  response  are  evident  at  <#»  =  ±5  deg.  Again, 
a  departure  from  quasisteady  behavior  is  seen  early  in  the  cycle 
(at  about  =  12  deg). 

When  the  static  (equilibrium)  flow  is  time- invariant,  the  NIR 
can  be  represented  as  the  difference  between  its  steady-state 
response  and  the  deficiency  function.^'^  The  deficiency  func¬ 
tion  is  defined  as 

F[<t>(^);  r,  T]  =  C;J<^(f);  oc.  r]  -  CJ0(f);  r,  r] 

where,  in  this  case,  the  steady-state  part  (first  term  on  the  RHS) 
is  the  derivative  of  the  static  rolling  moment  with  respect  to 
roll  angle.  (Also,  the  steady-state  term  on  the  RHS  can  be 
taken  as  the  mean  value  of  the  static  rolling-moment  derivative 
when  the  equilibrium  state  is  time  dependent.) 

Thus,  the  time- varying  part  of  the  NIR  is  represented  by  the 
deficiency  function.  When  this  form  of  the  NIR  is  put  into  the 
superposition  integral  [Eq.  (1)],  the  steady-state  part  integrates 
directly  to  give  the  change  in  the  static  coefficient  between  the 
limits  of  integration.  Therefore,  the  steady-state  parts  of  the 
initial  condition,  AC/,  and  the  two  integrals,  sum  to  the  instan¬ 
taneous  value  for  the  static  rolling-moment  coefficient 
C/.jt«[(#>(0].  Even  in  the  nonlinear  case,  with  critical  states  pres¬ 
ent,  the  total  response  under  dynamic  conditions  may  be  sep¬ 
arated  into  a  static  component  (evaluated  at  the  instantaneous 
roll  angle)  and  a  dynamic  component. 

The  dynamic  rolling  moment  coefficient  C/.^y,  is  shown  in 
Fig.  8d  for  a  roll-angle  offset  of  a  3-  and  12-deg  amplitude 
(the  same  motion  as  Fig.  8c).  Comparing  Figs.  8c  and  8d,  note 
that  the  static  and  dynamic  components  are  the  same  order  of 
magnitude.  Also,  as  shown  in  Fig.  8d,  locations  for  critical 
states,  identified  from  both  static  flowfield  structure  and  static 
force  behavior,  correlate  well  with  dramatic  changes  in  the 
dynamic  rolling-moment  response. 

The  results  discussed  earlier  for  a  limited  number  of  motions 
apply  over  a  wide  range  of  test  conditions  as  shown  in  Figs. 
9a-9c.  Each  of  these  is  a  contour  plot  of  the  dynamic  (total 
minus  static)  rolling-moment  coefficient  presented  in  the  phase 
plane.  Each  represents  k  =  0.08  data  for  the  complete  range  of 
test  amplitudes  at  a  given  roll-angle  offset.  A  series  of  tests 
with  fixed  offset  and  frequency  generates  a  family  of  ellipses, 
centered  about  the  offset  angle.  As  the  rolling  motion  pro¬ 
ceeds,  the  ellipses  are  traversed  in  the  clockwise  direction.  Off¬ 
sets  for  each  of  the  figures  are  0,  7,  and  14  deg,  respectively. 

The  contours  were  established  by  first  finding  a  mathemat¬ 
ical  representation  for  the  rolling-moment  data.  A  regression 
procedure,  similar  to  that  reported  in  Ref.  8,  was  used  to  do 
this.  A  correlation  coefficient  of  0.97  was  achieved  for  data 
covering  ail  test  conditions  at  Jt  =  0.08  and  0.14. 

In  all  cases  presented  in  Figs.  9a-9c,  the  contour  lines  turn 
rapidly,  becoming  essentially  parallel  to  the  roll-rate  axis  at 
0  =  ±5  deg.  Moving  clockwise  in  the  bottom  half  of  the 
ellipse  through  =  ~5  deg,  the  contour  lines  again  turn  rap¬ 
idly  between  -10  and  -15  deg  to  run  nearly  parallel  to  the 
roll  angle  axis.  Note  that  this  turning  point  is  less  distinct  than 
the  first,  perhaps  because  the  static  data  have  been  faired 
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Fig.  9  Dynamic  rolling-moment  contours:  a)  =  0,08,  4>o  =  0 
deg;  b)  ^  =  0.08,  =  7  deg;  and  c)  k  =  0.08,  4>o  =  14  deg. 

through  this  region.  The  pattern  repeats  in  the  upper  half  of 
the  figure.  Significantly,  the  loci  of  turning  points  are  indepen¬ 
dent  of  roll  rate,  further  evidence  that  they  are  critical  states 
because  they  are  determined  by  equilibrium  conditions. 

Following  a  critical-state  encounter,  the  dynamic  response 
becomes  weakly  dependent  on  roll  angle  (contour  lines  nearly 
parallel  to  the  roll- angle  axis).  Contours  in  this  region  are 
likely  the  result  of  expressing  time  variations  implicitly  in  the 
phase  plane  (rather  than  an  actual  dependence  on  roll  angle  or 
roll  rate). 

Finally,  note  that  as  the  offset  roll  angle  is  increased  (allow¬ 
ing  more  time  to  elapse  between  the  critical-state  encounter 
near  13  deg  with  positive  roll  rate  and  the  next  1 3-deg  critical- 
state  encounter  with  negative  roll  rate)  the  region  of  negligible 
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Fig.  10  Comparison  between  actual  data  and  locally  linear  model:  a)  A:  =  0.14,  =  3  deg,  and  =  5  deg  and  b)  A:  =  0.02,  =  3  deg, 

and  A<t>  ^  5  deg;  in*phase  and  quadrature  component  comparisons:  c)  =  0.14,  =  3  deg,  and  A<f>  =  5  deg  and  d)  A:  =  0.02,  i^o  =  3  deg, 

and  A<ff  =  5  deg;  e)  effect  of  extrapolating  low-frequency  derivatives  to  high  rates,  k  =  0.14,  <^»,  =  3  deg,  and  A<f>  =  18  deg;  and  f)  locally 
linear  model  vs  actual  data:  no  critical  states,  k  =  0.14,  =  28  deg,  and  =  5  deg. 


dynamic  rolling  moment  in  the  lower-right  quadrant  expands. 
There  is  more  time  for  the  transient  to  decay. 

Significance  of  Results 

Locally  linear  models  are  often  used  to  represent  aerody¬ 
namic  forces  and  moments  in  nonlinear  flight  simulations.  In 
this  model,  static  forces  and  moments  are  represented  by  non¬ 
linear  functions  of  instantaneous  values  of  angle  of  attack  and 
sideslip.  Dynamic  effects  are  calculated  by  using  locally  lin¬ 
ear  damping  derivatives  (linearized  about  the  instantaneous 


vehicle  state).  Dynamic  nonlinearities  are  accounted  for  by 
allowing  derivatives  to  be  functions  of  angle  of  attack  and 
sideslip. 

Linearization  of  the  damping  derivatives  is  often  a 
consequence  of  the  testing  technique  used  to  measure  them. 
Small-amplitude  harmonic  vehicle  motions  are  used  and  only 
the  quadrature  component  of  the  loads  at  the  frequency  of  the 
oscillation  is  considered.  Although  this  procedure  is  acceptable 
in  the  absence  of  significant  nonlinearities  (as  is  the  case  in 
many  small  amplitude  motions),  it  leads  to  results  that  cannot. 
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in  general,  be  legitimately  extrapolated  in  terms  of  rate,  which 
is  precisely  what  locally  linear  models  do. 

To  check  the  validity  of  locally  linear  damping  derivatives 
in  the  current  case,  forced  oscillation  data  were  reduced  in  this 
manner,  i.e.,  nonlinear  static  rolling-moment  data  were  used 
together  with  roll-damping  derivatives,  obtained  from  5-deg 
amplitude  tests  at  different  offset  angles.^®  Only  the  out-of- 
phase  rolling  moment  at  the  forcing  frequency  was  retained, 
although  up  to  20  harmonics  were  recorded.  Thus,  the  result 
is  a  derivative  linearized  about  the  mean  roll  angle  of  the  mo¬ 
tion.  The  locally  linear  model  was  then  used  to  predict  the 
measured  rolling  moment  (with  all  harmonics)  over  the  same 
motion  as  used  to  determine  the  damping  derivative. 

Comparisons  between  the  model  and  measured  responses 
are  shown  in  Figs.  10a- lOf.  In  Figs.  10a- lOd  and  lOf,  the 
damping  derivative  is  taken  to  be  constant  (consistent  with  the 
small  amplitude,  5  deg,  over  which  it  was  determined).  How¬ 
ever,  in  Fig.  lOe  the  damping  derivative  was  represented  by  a 
nonlinear  function  of  roll  angle  (owing  to  the  larger  amplitude, 
18  deg). 

Figures  10a  and  10b  show  results  for  <f>Q=  3  deg  and  reduced 
frequencies  of  it  =  0.14  and  0.02,  respectively.  The  correlation 
presented  in  Fig.  10a  is  totally  unacceptable.  Over  much  of 
the  cycle,  the  model  predicts  a  positive  rolling  moment,  while 
the  actual  response  is  the  opposite.  At  the  lowest  reduced  fre¬ 
quency  (k  =  0.02  in  Fig.  10b),  there  is  much  better  agreement, 
but  the  result  is  still  inadequate.  Even  at  this  extremely  low 
reduced  frequency,  there  is  a  lag  in  the  actual  response  that  is 
not  captured  by  the  locally-linear  model. 

Further  insight  is  provided  by  Figs.  10c  and  lOd,  which 
should  be  compared  to  Figs.  10a  and  10b,  respectively.  In- 
phase  (circular  symbols)  and  quadrature  (squares)  components 
of  the  actual  responses  are  compared  to  the  corresponding 
components  according  to  the  locdly-linear  model;  i.e.,  static 
rolling  moment  and  linearized  damping  terms,  respectively.  At 
it  =  0.14,  linearization  of  the  damping  term  (based  on  experi¬ 
mental  data  taken  at  this  frequency)  is  certainly  valid  over  the 
range  of  roll  rates  encountered  in  Fig.  10c.  However,  signifi¬ 
cant  errors  are  introduced  by  approximating  the  in-phase  com¬ 
ponent  with  the  static  data. 

Similarly,  at  k  =  0.02,  most  of  the  locally  linear  model’s 
error  is  again  seen  to  be  due  to  a  poor  prediction  of  the  in- 
phase  component  as  shown  in  Fig.  lOd.  However,  some  non¬ 
linear  effects  are  now  seen  in  the  quadrature  component.  The 
latter  problem  is  compounded  if  damping  derivatives  based  on 
small-amplitude,  low-frequency  (k  -  0.02)  data  are  used  to 
extrapolate  to  higher  angular  rates.  Such  a  case  is  shown  in 
Fig.  lOe  where  the  prediction  for  /:  =  0.14  with  3-deg  offset 
and  18-deg  amplitude  is  compared  with  the  corresponding  ob¬ 
served  rolling  moment.  Note  that  the  5-deg  amplitude,  the 
smallest  used  in  this  test  series,  is  two  to  three  times  larger 
than  that  used  in  typical  low-amplitude  tests.  Derivatives  based 
on  larger  amplitude  data  (with  their  associated  higher  rates) 
can  therefore  be  expected  to  yield  better  predictions  when  ex¬ 
trapolated  to  the  high  rates.  Even  so,  the  results  are  totally 
unacceptable;  the  in-phase  component  errors  are  completely 
swamped  by  those  in  the  damping  term.  Note  that  the  peaks 
in  the  prediction  are  centered  about  o)t  values  of  90  and  270 
deg  where  the  rate  is  the  highest 

Finally,  note  the  good  agreement  in  Fig.  lOf.  In  this  case, 
the  offset  roll  angle  of  28  deg  and  the  5-deg  amplitude  ensure 
that  the  critical  states  cannot  be  reached.  The  discrepancies 
here  are  likely  the  result  of  interpolation  errors  in  the  static 
component  Thus,  the  principal  difficulty  with  the  locally  linear 
model  noted  here  is  when  the  motion  includes  critical  states. 

Similar  effects  can  be  observed  in  Fig.  1 1  where  trajectories 
corresponding  to  three  free-to-roll  releases  are  superimposed 
on  a  phase-plane  portrait  constructed  with  the  locally  linear 
model.  Obviously  the  latter  is  totally  incapable  of  representing 
the  dramatic  motion  history  effects  and  nonlinearities  that  pre¬ 
vail  following  critical-state  encounters.  Therefore,  testing  tech- 


Fig.  11  Actual  free-to-roU  motions  vs  locally  linear  phase  por¬ 
trait 


niques  must  involve  model  motions  capable  of  eliciting  the 
transient  responses  not  observable  with  small-amplitude  low- 
rate  experiments.  Instantaneous  values  for  aerodynamic  loads 
must  be  measured  to  permit  observations  of  these  transients. 
Furthermore,  changes  in  skin-friction  and  flowfield  topology 
(together  with  a  careful  interpretation  of  static  data)  may  in¬ 
dicate  values  for  flow  and  motion  parameters  where  critical 
states  are  likely  to  occur. 

Clearly,  significant  errors  result  if  the  effects  of  critical-state 
encounters  are  not  handled  correctly.  Therefore,  aerodynamic 
models  that  can  account  for  the  existence  of  critical  states  and 
the  transients  induced  by  them  are  required. 
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Swept  and  delta  wings  maneuvering  at  moderate  and  high  angles  of  attack  produce  highly  nonlinear 
and  often  discontinuous  aerodynamic  forces  and  moments  that  are  difficult  to  model.  The  nonlinear 
indicia]  response  (NIR)  methodology  and  the  concept  of  critical  states  accompanied  by  changes  in  the 
flow  structure  and  topology  could  provide  a  rational  framework  for  the  analyses  and  modeling  of  these 
flows.  The  analysis  of  surface  oil-flow  photographs  and  laser  light  sheet  high-speed  video  images  of  smoke 
flow  has  been  performed.  The  correlation  of  the  structural  and  topological  changes  in  the  flow  with  force 
and  moment  data  follows.  Critical  states  are  often  accompanied  by  changes  in  the  flow  topology  and  not 
all  topological  changes  produce  measurable  changes  in  the  forces  and  moments,  however,  a  useful  rela¬ 
tionship  may  exist. 


Introduction 

IRCRAFT  with  delta  wings  have  flown  since  the  early 
1950s  (Vulcan,  F-102,  Mirage  IE,  MiG-21,  B-58,  XB-70, 
Eurofighter  2000  prototype,  etc.),  when  the  jet  engine  with  its 
higher  speed  capability  replaced  propellers  and  piston  engines. 
The  drag  reduction  for  thinner  wings  with  increased  sweep  was 
well  known  and  leads  to  the  delta  planfonn  with  greater  struc¬ 
tural  rigidity,  gradual  stall  characteristics,  and  high  lift  at  ma¬ 
neuvering  and  landing  attitudes  and  speeds.  The  triangular 
delta  wing  was,  however,  not  without  disadvantages  that  re¬ 
quired  stability  augmentation  and  double-delta,  straked  delta, 
and  other  closely  related  planforms  to  correct.  While  primarily 
used  for  high-performance  fighter  aircraft,  these  planforms  also 
have  found  application  to  the  Concorde  supersonic  transport 
and  proposed  High  Speed  Civil  Transport  designs. 

The  nonlinear  lift  increase  at  low  speeds  is  generated  by 
leeward  vortices  produced  by  the  rolling  up  of  the  shear  layers 
emanating  from  the  leading  edges.  These  vortices,  while  pro¬ 
ducing  the  desirable  lift  increase,  considerably  complicate  the 
flowfield  and  its  prediction.  The  flowfield  is  symmetric  and 
predictable  at  moderate  angles  of  attack  for  symmetric  flight 
conditions.  Unpredictable  changes  in  the  static  forces  and  mo¬ 
ments  occur  during  high  angle-of-attack  maneuvers  and  during 
even  moderate  asymmetric  flight  conditions  when  vortex  burst 
may  occur  over  the  planform.  Vortex  burst  or  vortex  break¬ 
down  is  characterized  by  the  sudden  expansion  of  the  highly 
organized  core  into  bubbles  or  spirals  along  the  core  axis. 
Shortly  downstream,  the  bubbles  or  spirals  diffuse  into  a  dis¬ 
organized,  swirling  turbulent  flow.  Vortex  burst  normally  oc¬ 
curs  first  in  the  wake,  proceeding  upstream  toward  the  wing’s 
apex  as  angle  of  attack,  sideslip  (windward  wing  only),  or 
aspect  ratio  increase.  Associated  with  the  forward  motion  of 
the  burst  point  are  loss  of  lift,  pitchup  and  nonlinear,  often 
discontinuous,  pitching  and  rolling  moment  characteristics.* 

Beginning  in  1987,  the  Canadian  Institute  for  Aerospace  Re¬ 
search  (lAR)  and  the  Flight  Dynamics  Directorate  of  Wright 
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Laboratory  have  been  pursuing  a  collaborative  research  pro¬ 
gram  on  unsteady  and  vortex-dominated  flows.  The  need  being 
addressed  is  the  lack  of  appropriate  mathematical  modeling 
techniques  to  represent  the  interaction  between  the  vehicle  mo¬ 
tion  and  the  forces  and  moments  created  by  these  flows.  The 
requirements  for  aerodynamic  modeling  arise  in  at  least  three 
areas:  1)  for  aerodynamic  understanding  and  design,  2)  for 
control  system  design,  and  3)  for  flight  simulation. 

A  theoretical  method  for  studying  the  nonlinear  aspects  of 
the  flight  dynamics  problem  has  been  under  development  by 
Tobak  and  his  colleagues^”'*  at  NASA  Ames  Research  Center 
since  the  early  1960s.  Their  approach  introduced  two  impor¬ 
tant  new  concepts:  1)  a  nonlinear  indicial  response  (NIR)  and 
2)  a  generalized  superposition  integral. 

The  NIR  approach,  which  has  been  derived  from  the 
Navier-Stokes  equations  for  a  time-invariant  equilibrium 
state,**  allows  that  critical  states  may  be  signaled  by  changes 
in  the  static  flow  topology,  often  manifest  in  the  position  and 
behavior  of  vortices  within  the  flow.  The  movement  of  the 
leading-edge  vortex  burst  point  onto  the  planform  and  the  in¬ 
troduction  of  secondary  and  tertiary  vortices  are  examples  of 
topological  changes,  as  are  changes  in  the  number  of  singular 
points  in  the  skin-friction  lines.  In  principle,  the  NIR  is  not 
restricted  to  a  time-invariant  equilibrium  state;  however,  the 
details  of  the  proof  have  not  been  completed  and  published. 

As  with  linear  indicial  response  methods,  the  arbitrary  mo¬ 
tion  is  represented  as  a  summation  of  responses  to  a  series  of 
step  motions.  The  NIR,  as  opposed  to  its  linear  counterpart, 
accounts  for  changes  induced  by  the  motion  history  leading 
up  to  step  onset.  Under  a  wide  variety  of  circumstances,  the 
summation  of  indicial  responses  approaches  the  generalized 
superposition  integral  in  the  limit. 

The  formulation  also  allows  for  critical  states  where  the 
aerodynamic  response  loses  analytical  dependence  on  the  mo¬ 
tion  variable,  such  as  when  aerc^ynamic  bifurcations  occur.^ 
A  critical  state  is  accommodated  by  splitting  the  generalized 
superposition  integral  at  the  critical  state  and  including  a  tran¬ 
sient  response  tenn.  Its  location  is  defined  by  the  values  of  a 
set  of  variables  that  characterize  the  instantaneous  motion. 

A  critical  state  is  a  flight  condition  where  large  and  persis¬ 
tent  transients  may  be  introduced  into  the  aerodynamic  loads, 
invalidating  the  linear  and  locally  linearized  aerodynamic 
models  traditionally  used  for  flight  mechanics  predictions.  For 
example,  the  result  of  roll-angle  and  roll-rate-induced  veloci¬ 
ties  at  the  leading  edge  of  a  delta  wing  is  an  unsymmetrical 
vortex  lift  that  produces  a  highly  nonlinear  (perhaps  discontin¬ 
uous)  change  in  the  rolling  moment. 

A  critical  state  is  often  associated  with  a  change  in  the  num¬ 
ber  of  critical  points  in  the  static  time-averaged  surface  flow 
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topology  as  shown  laten  Physically  realizable  flows  produce 
streamlines  and  skin  friction  line  patterns  in  oil  flows  that  must 
satisfy  topological  rules  based  entirely  on  kinematics  and  the 
differential  equations  of  streamlines.  A  stagnation  point  is  a 
critical  point,  also  called  singular  point  and  equilibrium  point/ 
in  the  flow  as  is  a  point  where  the  skin  friction  vector  is  zero. 
Critical  points  are  classified  as  nodes  or  saddles  depending  on 
their  appearance  and  mathematical  properties.  Nodes  are  fur¬ 
ther  subdivided  into  nodal  points  of  attachment  or  separation, 
foci,  and  centers.  Ail  skin-friction  lines  or  streamlines  intersect 
at  the  node  and  are  directed  away  from  a  node  of  attachment 
and  toward  a  node  of  separation.  The  focus  or  spiral  node  has 
all  of  the  lines  spiraling  inward  or  outward  around  the  critical 
point.  If  the  trajectories  form  closed  paths  around  the  focus,  it 
is  a  center.  The  saddle  has  two  oppositely  directed  lines  that 
intersect  at  the  critical  point.  All  other  lines  curve  away  in 
each  quadrant  to  avoid  this  critical  point. 

The  correspondence  between  critical  states  and  static  flow 
topology  changes  is  controversial,  especially  when  time-aver¬ 
aged  data  (e.g.,  surface  oil  flow  studies)  are  involved  (see  Ref. 
6).  The  analytic  dependence  of  the  aerodynamic  response  may 
be  lost  at  a  change  in  flow  topology,  i.e.,  when  a  change  in 
the  number  of  critical  points  is  observed.^  A  current  consensus 
is  that  critical  states  are  often  signaled  by  changes  in  the  flow 
topology;  however,  some  changes  in  the  flow  topology  do  not 
produce  measurable  force  and  moment  transients.  A  complete 
discussion  of  the  controversy  and  how  vortex  formation  in  the 
wake  of  a  circular  cylinder  at  low  Reynolds  number  provides 
a  counterexample  has  been  given  by  Jenkins  et  al.^ 

Experimental  Program 

The  flow  about  delta  wings  is  highly  sensitive  to  very  small 
differences  in  model  geometiy,  motion  variables,  and  flow  con¬ 
ditions  making  analysis  and  correlation  of  the  vast  amounts  of 
available  data  difficult.  The  experimental  program  eliminated 
many  of  these  variables  by  testing  the  same  model  in  two  wind 
tunnels  at  similar  flow  conditions  with  overlapping  ranges  of 
the  motion  variables  and  different  support  systems.  The  pro¬ 
gram  is  continuing  with  previous  measurements  identifying 
regions  of  particular  interest  for  succeeding  tests. 

The  common  model  is  a  65-deg  swept  triangular  wing  with 
sharp  leading  edges  (10-deg  symmetric  nonconical  bevel  on 
all  three  edges)  and  centerbody  (Fig.  1).  It  was  especially  de¬ 
signed  and  constructed  to  be  lightweight,  with  low  inertia  to 
attain  the  desired  motions,  yet  strong  enough  to  withstand  the 
loads  encountered  at  the  high  pitch  and  roll  rates  required  to 
match  full-scale  aircraft  reduced  frequencies  (see  Hanff  and 
Jenkins’  for  details). 

Static  and  dynamic  tests  of  this  model  have  been  a  partic¬ 
ularly  rich  source  of  nonlinearities  and  other  unusual  dynamic 
behavior.  The  extensive  experimental  database  allows  concur¬ 
rent  study  of  vortex  dynamics  and  the  resulting  unsteady  aero- 

Composite  Construction 


(0.58  m) 

Fig.  1  Model  geometry. 


dynamic  forces  and  moments.  The  correct  identification  of 
changes  in  the  flow  physics  and  topology  is  crucial  to  further 
understanding  the  NIR  methodology  and  its  application  to  this 
series  of  wind-tunnel  test  data. 

The  lAR  experimental  program,  conducted  in  the  National 
Aeronautical  Establishment  6  X  9  ft  low-speed  wind  tunnel, 
has  been  described  by  Hanff  and  Jenkins’  and  Huang  et  al.‘ 
The  Wright  Laboratory  7  X  10  ft  Subsonic  Aerodynamic  Re¬ 
search  Laboratory  and  the  experimental  data  has  been  de¬ 
scribed  by  Jenkins  and  various  coauthors.**®’’  LAR  personnel 
conducted  all  dynamic  testing  using  their  roll  rig.  The  most 
extensive  data  set  was  taken  at  Mach  number  0.3,  Reynolds 
number  3.6  X  10*  based  on  the  centerline  chord  of  2.04,  (0.622 
m)  sting  angle  of  30  deg,  with  various  static  and  mean  roll 
angles  <f>y  and  roll  rates.  The  static  roll  response  is  highly  non¬ 
linear, with  three  distinct  static  roll  attractors.  Further,  forced- 
oscillation  and  free-to-roll  data  exhibited  strong  nonlinearities 
that  cannot  be  modeled  using  locally  linear  models.* 

Critical  States  and  Topology 

Analysis  of  the  static  and  dynamic  force  and  moment  data* 
indicated  critical  states  near  <^  =  5  and  11.3  deg,  based  on 
slope  changes  and  discontinuities  in  the  static  time-averaged 
pitching  and  rolling  moments  and  transients  in  the  dynamic 
data.  Additional  evidence  confirming  the  existence  of  these 
critical  states  has  been  provided  by  a  wealth  of  subsequent 
experimental  data®*’  and  analyses.**’^  The  following  analysis 
correlates  the  changes  in  static  forces  and  moments  with  the 
topology  changes  described  by  Huang  et  al.*  The  changes  in 
the  forces  are  small  because  the  vortex  flow  patterns  contribute 
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Fig.  3  Pitching  moment  coefficient:  a)  full  range,  —  90  ^  ^  90 

deg  and  b)  expanded  scale,  0  ^  20  deg. 

less  than  10%  of  the  total  lift  at  30-deg  angle  of  attack,  zero 
roll  angle  for  this  aspect  ratio  1.86  wing.  Thus,  changes  in  the 
vortex  flows  rarely  obviate  trends  established  by  the  attached 
flow,  as  will  be  discussed  subsequently.  Criticd  states,  how¬ 
ever,  lead  to  large  deviations  from  the  static  characteristics, 
especially  in  the  moments. 

Experimental  data  for  the  entire  roll  angle  range  (-~90  to 
+90  deg)  and  to  an  expanded  scale  emphasizing  the  5  and 
11.3  deg  (to  be  discussed  later)  critical  states  are  shown  in 
Figs.  2-5.  The  force  and  moment  coefficients  are  based  on 
the  dynamic  pressure  and  the  wing  area  of  1.94  ft^  (0.180  m^). 
Reference  lengths  are  the  mean  aerodynamic  chord  of  1.36  ft 
(0.415  m)  for  the  pitching  moment  coefficient  and  the  wing 
span  of  1.90  ft  (0.58  m)  for  the  body-axis  rolling  moment 
coefficient.  Dashed  vertical  lines  indicate  the  critical  state  lo¬ 
cations.  The  data  has  been  sorted  and  plotted  according  to  the 
direction  the  wing  was  rolling  toward  (increasing  or  decreasing 
cfi)  when  the  particular  roll  angle  was  attained.  This  is  contrary 
to  previous  analyses*®’’  based  on  rather  sparse  static  data  and 
permits  a  more  detailed  examination  for  discontinuities,  slope 
changes,  and  hysteresis  that  indicate  the  crossing  of  critical 
states. 

The  full-range  data  plots.  Figs.  2a,  3a,  4a,  and  5a,  show 
proper  symmetry.  Time-averaged  force  and  moment  behavior 
at  the  previously  determined  critical  states  at  large  roll  angles’ 
(-51.3,  -49.5,  50.1,  and  51.4  deg),  are  also  shown  on  these 
figures  as  dashed  vertical  lines.  The  data  density  on  this  scale 
precludes  definite  identification  of  the  critical  states  at  either 
<^  =  5  or  11.3  deg.  All  data  were  fit  using  Legendre  polyno¬ 
mials  and  stepwise  regression  analysis  as  reported  previously*’’ 
on  the  earlier  sparse  data  sets. 


Five-Deg  Critical  State 

The  discontinuity  near  5-deg  roll  angle  in  Figs.  2b,  3b,  and 
4b  is  clearly  a  critical  state.  It  is  associated  with  the  leeward 
(port  for  positive  roll  angles)  wing  vortex  burst  rapidly 
crossing  the  trailing  edge  during  extremely  small  roll  angle 
increases,  as  shown  by  previous  authors.*’®’’  This  is  in  agree¬ 
ment  with  the  topologicd  maps  of  Huang  et  al.‘  that  show  the 
primary  vortex  burst  is  over  the  leeward  wing  at  =  4  deg 
and  is  absent  at  <^  =  7  deg.  It  is  only  the  secondary  and  tertiary 
vortices  that  lift  off  [shown  as  whorls  or  foci  (spiral  nodes)  in 
the  skin-friction  topology]  at  <^  =  7  deg. 

Figure  6a  shows  the  vortex  burst  of  the  leeward  leading- 
edge  vortex  is  10-15%  chord  aft  of  the  trailing  edge  at  ^  = 
—5.4  deg,  while  Fig.  6b  shows  the  vortex  burst  has  moved 
upstream  onto  the  planform  at  4>  =  -5.0  deg.  Note  the  rapid 
forward  progression  to  10-15%  chord  ahead  of  the  trailing 
edge  in  only  0.4-deg  roll  angle  change.  This  is  the  discrete 
change  that  defines  the  critical  state. 

The  increase  in  positive  rolling  moment  as  <f)  increases  (Fig. 
2b)  is  consistent  with  the  formation  of  the  coherent  vortex  over 
the  entire  leeward  wing.  As  the  vortex  burst  point  moves  aft, 
off  the  wing,  the  extended,  concentrated  vortex  core  decreases 
the  pressures  on  the  leeward  wing  upper  surface  and  contrib¬ 
utes  an  additional  destabilizing  moment.  The  additional  force 
near  the  trailing  edge  also  contributes  a  nose-down,  stabilizing 
increment  to  the  pitching  moment  (Fig.  3b).  The  axial  force 
(Fig.  4b)  shows  an  increase,  consistent  with  the  known  drag 
increase  attendant  with  lifting  leading-edge  vortices.  The  nor¬ 
mal  force  is  decreasing  because  of  the  sideslip  increase. 

The  regression  analysis  performed  on  data  in  this  region  was 
split  into  two  separate  domains,  increasing  and  decreasing  </>, 


Fig.  4  Axial  force  coefficient:  a)  full  range,  —  90  ^  90  deg 

and  b)  expanded  scale,  0  ^  <l>  20  deg. 
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Fig.  5  Normal  force  coefficient:  a)  full  range,  —90  <  <  90  deg 

and  b)  expanded  scale,  0  ^  20  deg. 


between  0  =  5. 0-5. 7  deg.  Only  the  increasing  <p  data  were 
used  in  this  region  for  the  regression  fitting  the  lower  branch 
in  this  region,  while  only  decreasing  (f>  data  were  used  for  the 
upper  branch.  This  partitioning  of  data  is  consistent  with  the 
definition  of  a  static  hysteresis,  and  shows  that  a  narrow  hys¬ 
teresis  loop  apparently  exists  in  this  region.  Further,  two  crit¬ 
ical  states  must  exist  here,  one  at  the  terminus  of  each  branch. 
Past  and  present  analyses  of  experimental  data  suggest  such  a 
possibility.  The  analysis  of  Addington  and  Jenkins®  show  that 
the  leeward  leading-edge  vortex  burst  point  crosses  the  trailing 
edge  from  being  over  the  planform  at  approximately  <p  =  5 
deg.  Huang  et  al.*  analysis  confirms  this  finding,  as  well  as 
suggests  that  an  additional  topological  change  occurs  in  this 
region. 

Several  correlations  have  been  devised  to  predict  the  move¬ 
ment  of  vortex  burst  points  with  leading-edge  sweep  angle, 
angle  of  attack,  and  roll  angle.  The  predicted  roll  angle  for 
vortex  breakdown  at  the  trailing  edge  is  7  deg  based  on  a  linear 
regression  fit®  of  the  static  data  of  Hanff  and  Huang^*;  8.8  deg 
based  on  Ericsson  and  Hanff’ s'®  method  (using  Wentz  and 
Kohlman’s^^  vortex  breakdown  correlations,  vortex  breakdown 
occurs  at  the  trailing  edge  at  19-deg  angle  of  attack  for  <^  =  0 
deg);  and  13.7  deg  from  the  revised  Huang  and  Hanff^** 
method.  Their  correlations  are  a  better  approximation  to  ex¬ 
perimental  data'^’'*^  for  the  60-  and  70-deg  swept  wings  than 
for  the  65-deg  swept  wing,  again  indicating  unusual  behav- 
ior^'‘°  at  this  combination  of  roll  and  sweep  angles.  None  of 
these  vortex  breakdown  prediction  methods  account  for  the 
rapid  jump  in  vortex  position  because  they  assume  smooth 
functions  to  correlate  the  experimental  data. 


11.3-Deg  CriticaJ  State 

The  tentatively  associated  change  in  flow  topology  near  (t> 
=  12  deg,  the  vortex  breakdown  point  reaching  the  windward 
(starboard)  wing’s  apex,®  was  shown  to  be  erroneous®  by  anal¬ 
ysis  of  the  subsequent  static  and  forced-oscillation  dynamic 
wind-tunnel  data  called  for  by  Jenkins  et  al.®  The  assumption 
of  vortex  breakdown  reaching  the  apex  of  the  windward  wing 
near  (p  -  12  deg  was  not  unfounded,  particularly  in  light  of 
the  lack  of  detailed  flow  visualization  data  available  at  the 
time.  The  linear  regression  data  fit  of  Addington  and  Jenkins® 


Fig.  6  Leeward  leading-edge  vortex  burst:  a)  in  the  wake,  </>  = 
-5.4  deg  and  b)  over  the  planform,  =  —5.0  deg. 


Fig.  7  Laser  light  sheet  flow  visualization. 
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Fig.  8  Surface  oil  flow,  — 7  deg. 


predicted  =  13  deg,  whereas  the  methods  of  Ericsson  and 
Hanff*°  and  Huang  and  Hanff^'*  predicted  12.6  and  12.15  deg, 
respectively.  Vortex  breakdown  would  occur  aft  of  the  trailing 
edge  on  the  leeward  wing  and  could  not  cause  the  observed 
force  and  moment  transients.  The  preceding  analysis  methods 
generally  agree  with  the  original  association  of  vortex  break¬ 
down  at  the  apex.  However,  analysis  of  the  previous  dynamic 
and  new  high-resolution  static  data  taken  in  the  same  wind 
tunnel  using  the  same  model  (the  sting  and  instrumentation 
were  different,  which  permitted  additional  verification  of  the 
data  through  correlations  at  selected  tie-in  conditions)  clearly 
shows  that  the  critical  state  is  not  associated  with  vortex  break¬ 
down  advancing  to  the  apex  at  <f>  =  11.3  deg. 

The  experimental  setup  with  the  laser  light  sheet  is  shown 
in  side  view  in  Fig.  7.  The  laser  light  sheet  has  progressed  to 
aft  of  the  midchord  and  the  leading-edge  vortex  over  the  entire 
leeward  wing  is  evident  at  </>  =  14.5  deg.  The  light  sheet  shows 
the  shear  layer  from  the  windward  wing  has  moved  across  the 
centerbody  at  this  chordwise  station,  in  agreement  with  the 
observations  of  Hsia  et  al.^  \fideo  recordings  of  the  fore  and 
aft  sweeps  of  the  laser  light  sheets  show  that  the  flow  becomes 
unsteady  as  the  critical  roll  angle  is  approached  and  returns  to. 
steady  motion  as  this  roll  angle  is  exceeded.  While  a  still  im¬ 
age  suitable  for  publication  could  not  be  captured  by  our 
equipment,  the  videotapes  clearly  show  a  vortex  near  the  apex 
of  both  the  windward  and  leeward  wings  at  <^>  =  15  deg.  Both 
vortices  are  evident  at  all  roll  angles  less  than  15  deg,  while 
the  vortex  over  the  windward  wing  is  not  apparent  ai  (f>  =  18 
deg.  The  vortex  persists  over  the  entire  leeward  wing  at  this 
roll  angle. 

Surface  oil  flows  near  the  trailing  edge  at  nearly  the  same 
test  conditions  (30-deg  roll  angle  and  0.3  Mach  number)  are 
shown  in  Figs.  8  and  9  at  roll  angles  of  -7  and  — 14  deg.  A 


description  of  the  model  and  experimental  setup  is  contained 
in  Hanff  et  al.*^  The  change  in  flow  topology  is  evident.  Two 
vortex  whorl  patterns  are  evident  over  the  leeward  wing  at  a 
roll  angle  of  —7  deg,  while  one  is  evident  at  =  — 14  deg. 
These  patterns  correspond  to  the  liftoff  of  the  tertiary  vortex 
(Fig.  8)  and  to  liftoff  of  both  the  secondary  and  tertiary  vor¬ 
tices  (Fig.  9)  on  the  planform.  The  flow  pattern  and  topology 
sketches  from  Huang  et  al.^  show  this  topology  change  occurs 
between  roll  angles  of  ~7  and  — 14  deg,  but  the  precise  roll 
angle  and  triggering  mechanism  are  unknown.  According  to 
Huang  et  al.,‘  the  change  in  flow  pattern  corresponds  to  a 
critical  state  crossing;  however,  additional  images  and  analysis 
are  needed  to  fill  these  experimental  data  gaps. 

From  Huang  et  al.,^  between  10-  and  14-deg  roll  angle,  the 
flow  on  the  windward  wing  changes  from  an  organized  vortex 
flow  with  the  tertiary  vortex  lifting  off  the  surface  to  a  topol¬ 
ogy  representative  of  primary  vortex  breakdown  advancing  to¬ 
ward,  but  not  attaining,  the  wing  apex.  This  creates  a  disor¬ 
ganized  spiraling  flow  aft  of  the  breicdown,  exhibiting  reverse 
flow  partly  swept  across  from  the  leeward  wing,  and  a  loss  of 
vortex  lift  on  the  windward  wing.  The  leeward  wing  imdergoes 
a  topological  change  indicative  of  reattachment  of  the  second¬ 
ary  vortex  near  the  wingtip  as  the  roll  angle  increases.  The 
organized,  spiraling  flow  beneath  this  secondary  vortex  in¬ 
creases  the  normal  force  on  the  area  near  the  wingtip  and  pro¬ 
duces  force  and  moment  changes  similar  to  the  5-deg  case, 
albeit  smaller.  The  force  and  moment  changes  are  also  influ¬ 
enced  by  the  vortex  lift  loss  over  the  windward  wing. 

Returning  to  Figs.  2-5,  at  =  1 1.3  deg,  the  negative  slope 
of  the  rolling  moment  curve  becomes  less  negative  because  of 
the  additional  destabilizing  moments  contributed  by  both  the 
windward  and  leeward  wings  (Fig,  2b).  The  nose-up  pitching 
moment  slope  is  increased  by  the  additional  force  near  the 
leeward  wing  trailing  edge  (Fig.  3b),  causing  the  discontinuity 
in  slope.  The  scatter  in  the  axial  force  data  (Fig.  4b)  and  the 
magnitude  of  the  changes  precludes  definitive  conclusions; 


Fig.  9  Surface  oil  flow,  =  ”14  deg. 
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however,  a  discontinuity  in  slope  is  evident.  The  overall  nor¬ 
mal  force  is  decreasing  (Fig.  5b)  because  of  increasing  roll 
angle  and  the  additional  force  caused  by  the  leeward  vortex  is 
somewhat  counteracted  by  the  reduced  vortex  lift  over  the 
windward  wing.  The  rate  of  lift  decrease  is  slowed  and  causes 
the  slope  discontinuity. 

8.5-Deg  Critical  State 

A  suspected  critical  state  is  marked  at  <f>-  8.5  deg  in  Figs. 
2b-5b.  The  topological  change  shown  by  Huang  et  al.^  rein¬ 
forces  the  evidence  obtained  by  separately  fitting  the  roll  angle 
increasing  or  decreasing  data. 

The  polynomial  regression  fit  of  the  experimental  force  and 
moment  data  for  increasing  and  decreasing  roll  angle  is  shown 
as  solid  lines  in  Figs.  2b,  3b,  4b,  and  5b.  Two  distinct  branches 
near  =  8.5  deg  are  evident  in  each  data  set.  The  use  of  an 
overlapping  domain  (from  (j>  =  7.7  to  9.2  deg)  to  fit  the  data 
increases  the  correlation  coefficient  compared  to  strictly  mono¬ 
tonic  data  fits. 

Huang  et  al.‘  report  that  a  topology  change  occurs  between 
<t>-l  and  10  deg.  This  change  is  associated  with  global  flow 
separation  occurring  over  the  windward  (starboard)  wing  as 
roll  angle  increases.  The  force  and  moment  discontinuities  oc¬ 
curring;  increase  in  rolling  moment,  decrease  in  nose-up  pitch¬ 
ing  moment,  and  decrease  in  axial  and  normal  forces  are  con¬ 
sistent  with  the  loss  of  suction  because  of  the  disappearance 
of  the  leading-edge  vortex  over  the  windward  wing. 

Discussion 

The  previous  correlations  agree  with  the  topological  analysis 
of  Huang  et  al.^  that  is  based  on  oil  flow  pictures  and  video 
images  taken  at  discrete  roll  angles  every  3  or  4  deg.  Addi¬ 
tional  anomalies  in  the  force  and  moment  data  create  the  pos¬ 
sibility  of  additional  topological  changes.  The  separate  analysis 
of  roll  angle  increasing  or  decreasing  data  indicates  the  pos¬ 
sibility  of  hysteresis  at  the  5-deg  critical  state  (Fig.  2b)  that 
would  require  two  critical  states,  one  near  5,0  deg  and  one 
near  5.7  deg.  Hysteresis  is  also  evident  in  the  pitching  moment 
data  (Fig.  3b)  and  an  additional  discontinuity  appears  near  <f> 
=  8  deg,  where  flow  visualization  data  were  not  recorded. 

A  presently  unexplained  anomaly  appears  in  these  data.  The 
branches  in  the  force  and  moment  data  at  <^  =  .8.5  deg  cannot 
form  a  hysteresis  loop.  The  increases  or  decreases  in  the  co¬ 
efficients  are  either  in  the  opposite  direction  or  at  opposite  ends 
of  where  they  would  occur  in  a  hysteresis  loop.  Furthermore, 
only  one  topology  change  is  observed  between  <^  =  7-10  deg 
in  the  mean  surface  oil  flows.  The  lack  of  data  in  this  roll 
angle  range  does  not  permit  resolution  of  this  anomaly,  how¬ 
ever,  several  possible  explanations  are  1)  additional  topologi¬ 
cal  changes  occur  in  the  surface  flow,  2)  topological  changes 
occur  in  the  streamlines  of  the  off-body  flow,  or  3)  the  flow 
becomes  unsteady  through  a  Hopf  bifurcation. 

Current  analyses  are  based  on  surface  oil  flows  related  to 
static  force  and  moment  data.  The  time-averaging  property  of 
this  flow  visualization  method  precludes  inference  of  the  un¬ 
steady  effects  necessary  to  farther  elucidate  the  relationship 
between  critical  states  in  the  aerodynamic  response  and 
changes  in  the  number  of  critical  points  in  the  flow. 

Certainly  a  wealth  of  interesting  possibilities  occurs  that 
cannot  be  confirmed  based  on  presently 'available  data,  but 
point  to  the  need  for  specific  high  resolution  data  to  resolve 
these  questions.  The  high  resolution  available  from  Navier- 


Stokes  numerical  simulations*®  of  these  experiments  cannot  be 
equaled  for  directing  attention  to,  and  clarifying,  otherwise 
confusing  or  unnoticed  flow  features.  It  is  currently  too  ex¬ 
pensive  to  perform  simulations  at  all  of  the  points  needed.  The 
need  for  additional  testing  to  fill  these  gaps  is  evident. 

Correlation  of  critical  states  and  topology  changes  is  ex¬ 
ceedingly  difficult  and  time  consuming  because  of  the  preci¬ 
sion  and  resolution  required  in  both  the  static  and  dynamic 
data.  The  static  force  and  moment  data  across  suspected  dis¬ 
continuities  and  hysteresis  loops  should  be  used,  along  with 
flow  visualization,  to  identify  conditions  for  critical  states. 
Subsequent  dynamic  experiments  should  focus  on  these  con¬ 
ditions  to  identify  critical-state  transient  characteristics. 
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Introduction 

The  emergence  and  consequences  of  asymmetries  in  swirling 
flows  that  are  initially  steady  and  axisymmetric  are  exam¬ 
ined.  The  strength  of  an  isolated  vortex  in  a  tube  is  increased  in 
a  parametric  fashion  through  a  critical  value,  where  stability  to 
three-dimensional  disturbances  is  lost.  The  flow  behavior  under¬ 
goes  a  bifurcation  at  the  critical  value  from  steady  and  axisymmetric 
flow  to  unsteady  and  three-dimensional  flow.  Other  computations 
of  bifurcation  phenomena  in  swirling  flows  have  been  presented 
by  Leibovich  and  Kribus,^  Beran  and  Culick,^  and  Lopez.^  These 
works  are  limited  to  bifurcations  that  only  involve  axisymmetric 
flows. 

Axisymmetric  base  flows  serve  as  initial  conditions  to  a  three- 
dimensional  time-integration  algorithm.  The  minimum  axial  veloc¬ 
ity  component  Q(0  is  computed  and  compared  with  the  initial  value. 
Of  particular  interest  is  the  characterization  of  the  stability  loss  and 
the  relationship  between  the  appearance  of  asymmetries  and  the 
associated  changes  in  Q, 

The  computational  approach  is  as  follows.  First,  a  pseudo- 
arclength  continuation  (PAG)  algorithm^  provides  the  steady,  ax¬ 
isymmetric  initial  condition  for  a  specified  vortex  strength  V.  The 
Mach  number  M  and  Reynolds  number  Re  (based  on  vortex  core  ra¬ 
dius)  are  held  fixed  at  0.3  and  2.5  x  10^,  respectively.  No  nonunique 
axisymmetric  solutions  are  found  at  Re  =  2.5  x  10^,  consistent 
with  Ref.  2,  The  two-dimensional  solution  is  then  interpolated 
onto  the  three-dimensional  mesh  using  a  fourth-order-accurate  cu¬ 
bic  spline  scheme,'*  Then  time  integration  is  carried  out  by  the 
time-accurate  Navier-Stokes  (TANS)  model.  The  TANS  model  is  a 
special-purpose,  time-integration  algorithm  developed  specifically 
for  this  work  and  is  described  in  Ref.  5.  The  TANS  model  em¬ 
ploys  fourth-order  compact,  or  Fade,  operators*'  to  discretize  spa¬ 
tial  derivatives,  thus  allowing  for  fewer  grid  nodes  while  main¬ 
taining  sufficient  accuracy.  A  multiblock  grid  is  used  to  allow  for 
a  nearly  rectilinear  arrangement  of  nodes  near  the  tube  center- 
line,  while  near  orthogonality  is  maintained  at  the  tube  wall.  The 
PAG  algorithm  is  implemented  with  the  same  boundary  conditions 
and  tube  geometry  as  the  TANS  model,  using  a  simple  algebraic 
grid. 

The  physical  domain  consists  of  a  two-stage  cylindrical  tube  of 
circular  cross  section  and  varying  radius.^  The  first  stage  contains  a 
constriction  that  controls  the  upstream  movement  of  the  breakdown 
region.  The  tube  radius  (nondimensionalized  by  vortex  core  radius) 
at  the  inlet  station  is  fixed  at  Ro  =  2.  The  number  of  nodes  in 
the  computational  coordinate  directions  are  (nx,  ny,  nz),  where  nx 
defines  streamwise  spacing  and  ny  and  nz  are  equal  and  define 
cross-plane  spacing  in  the  y  and  z  directions,  respectively.  Three 
grids  are  employed  in  this  work.  Grid  G1  consists  of  98  x  41^ 
nodes,  grid  G2  contains  122  x  61^  nodes,  and  gridG3  uses  146x41^ 
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nodes.  The  nondimensional  tube  length  L  for  grids  G1  and  G2  is 
20,  whereas  the  tube  length  for  grid  G3  is  30.  Grid  G2  uses  axial 
clustering  in  the  breakdown  region  to  achieve  a  minimum  axial 
spacing  of  about  0. 1  at  x  =  5.  The  time  step  At  for  runs  using  grids 
GI  and  G3  is  0.04,  whereas  runs  on  the  finer  grid  G2  have  a  time 
step  of  0.025. 

Axisymmetric,  columnar  conditions  are  specified  at  the  inflow 
boundary  plane.  The  inlet  axial  velocity  is  chosen  to  be  uniform. 
The  swirl  velocity  component  v  is  modeled  from  a  Burger  vortex 
and  is  appropriate  for  modeling  the  profiles  obtained  from  a  swirl- 
vane  apparatus’; 

C(0,r.e)  =  Vr-' (1) 

where  r  and  9  denote  the  radial  and  azimuthal  coordinate  directions 
and  V  is  the  vortex  strength.  A  nonuniform  inflow  density  profile 
based  on  columnar  flow  is  prescribed  as  a  Dirichlet  condition.^  Slip 
is  allowed  along  the  tube  wall  surface,  yielding  an  impermeabil¬ 
ity  boundary  condition  that  relates  the  axial  and  radial  velocities. 
Outflow  conditions  reflect  an  assumed  columnar  flow  state. 


Results 

Between  V=  1.0  and  1.5,  the  initially  steady,  two-dimensional 
flow  does  not  evolve  in  time  to  an  asymmetric  flow  state.  However, 
as  V  is  increased  slightly  to  1 .53,  the  flow  becomes  three  dimensional 
and  time  periodic.  As  V  is  further  increased,  the  magnitude  of  the 
asymmetry  also  increases.  This  change  in  stability  is  attributed  to 
the  crossing  of  a  supercritical  Hopf  bifurcation  point  somewhere 
between  V  =  1.5  and  1.53. 

The  deviations  from  the  two-dimensional  initial  conditions  are 
illustrated  in  Fig.  1,  where  Q  (solid  lines  of  Figs,  la- Id)  is  plotted 
vs  time  for  V  =  1.5,  1.53,  1.55,  and  1.65.  In  Fig.  la  (V  =  1.5),  it 
is  evident  that  no  appreciable  deviation  from  the  initial  condition  is 
present.  This  is  confinned  by  monitoring  the  maximum  solution  cor¬ 
rection  at  each  iteration.  At  V  =  1 .53  (Fig.  lb),  the  solution  departs 
from  the  initial  condition,  as  evident  by  the  increase  in  Q  as  time 
increases.  The  term  Q(t)  eventually  becomes  time  periodic,  but  this 
behavior  is  not  discernible  from  the  scales  of  the  figure.  Periodicity 
is  confirmed  using  phase  plots  of  the  three  velocity  components  at 
a  fixed  centerline  node  near  x  =  7.  At  V  =  1 .55  and  1.65  (Figs.  Ic 
and  Id),  much  larger  deviations  are  evident. 

The  effect  of  grid  resolution,  time  step,  and  tube  length  are  also 
illustrated  in  Fig.  1.  Figures  la-lc  show  how  grid  refinement  from 
Gl  to  G2  affects  Q.  At  V  =  1.5,  resolving  the  grid  from  Gl  to 
G2  results  in  a  slight  decrease  in  the  time-asymptotic  value  of  Q. 
However,  the  flow  remains  steady  and  axisymmetric.  At  V  =  1.53, 
grid  refinement  results  in  a  more  pronounced  transient  response. 
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Fig.  1  Temporal  behavior  of  Qt  grid  and  time-step  sensitivity  and  tube 
length  L  sensitivity:  a)  V  =  1.5  (Gl:  coarse  grid,  G2:  fine  grid,  L  -  20), 
b)  V  =  1.53,  c)  V  =  1.55,  d)  V  =  1.65,  e)  effect  of  time  step  for  V  =  1.55, 
and  f )  effect  of  tube  length  for  V  =  1.53  (G3:  coarse  grid,  L  =  30). 
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possibly  representing  a  slight  shift  of  the  Hopf  point  to  smaller 
values  of  vortex  strength.  The  sensitivity  to  grid  refinement  appears 
to  diminish  as  vortex  strength  is  increased  away  from  the  Hopf 
point,  as  evident  in  Fig.  Ic.  Figure  le  shows  that  reducing  the  time 
step  from  0.04  to  0.025,  with  grid  Gl,  has  a  negligible  influence  on 
the  transient  behavior  of  Q.  This  suggests  that  the  differences  in  Q 
evident  in  Figs,  la-lc  are  due  primarily  to  grid  refinement  and  not 
to  the  reduction  of  the  time  step.  Finally,  the  effect  of  tube  length  is 
shown  in  Fig.  If,  where  increasing  the  tube  length  I  from  20  to  30 
appears  to  have  a  small  but  negligible  effect  on  Q. 

The  nature  of  the  Hopf  bifurcation  is  illustrated  in  Fig.  2.  To  iden¬ 
tify  the  onset  of  three-dimensionality,  a  global  parameter  H  is  con¬ 
structed  and  is  defined  to  be  the  maximum  absolute  value  of  dv/d^. 
The  term  H  is  nonzero  when  the  flow  is  asymmetric.  In  Fig.  2a, 
iY  abruptly  departs  from  zero  between  V  =  1.5  and  1.53,  Within 
this  range  of  V,  fiilly  three-dimensional  solutions  bifurcate  from 
the  branch  of  two-dimensional  solutions  when  the  two-dimensional 
solutions  become  physically  unstable.  Flow  unsteadiness  and  asym¬ 
metry  are  characterized  by  another  parameter  Q,  which  is  defined 
as  the  minimum  and  maximum  values  of  the  sum  of  the  cross-plane 
velocity  components  along  the  tube  centerline.  By  definition,  Q  is 
zero  for  an  axisymmetric  flow.  In  addition,  since  the  cross-plane  ve¬ 
locity  components  are  unsteady,  periodic  functions  (past  V  =  1.5), 


g.  3  Particle  traces  showing  development  of  three-dimensionality  with  increasing  vortex  strength:  a)  V  =  1.5,  b)  V  =  1.53,  c)  V  =  1.65,  and  d) 
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Q  characterizes  the  degree  of  flow  unsteadiness  by  representing  the 
largest  amplitude  of  the  cross-plane  velocity  sum  along  the  center¬ 
line.  Rgure  2b  shows  that  Q  departs  from  zero  at  the  same  value  of 
V  that  the  flow  becomes  asymmetric.  Thus,  Figs.  2a  and  2b  demon¬ 
strate  that  flow  unsteadiness  and  asymmetry  are  intimately  linked. 
The  loss  of  stability  to  time-periochc  flow  is  evidence  for  a  Hopf 
bifurcation.  The  bifurcation  is  supercritical,*  since  the  amplitude 
of  the  disturbance,  characterized  by  Q,  grows  from  zero  as  V  is 
increased  past  the  Hopf  point.  The  initial  and  final  values  of  Q 
are  shown  in  Fig.  2c.  A  region  in  V  exists  where  two-dimensional 
flows  with  bubble-type  breakdown  evolve  into  three-dimensional, 
unsteady  flows  with  no  reversed  flow.  Also,  it  is  of  interest  to  note 
that  the  Hopf  point  and  the  appearance  of  reversed  flow  occur  at 
different  values  of  V.  Thus,  loss  of  stability  is  not  a  consequence  of 
a  gross  structural  change  in  the  flowfield.  Other  examples  in  which 
changes  in  flowfield  topology  are  disassociated  from  hydrodynamic 
bifurcations  are  given  by  Lopez.^ 

The  observed  increase  in  Q  as  the  flow  transitions  from  two- 
to  three-dimensional  flow  is  correlated  to  a  particular  asymmet¬ 
ric  term  in  the  governing  equations.  Details  of  this  analysis  can 
be  found  in  Ref.  5.  The  correlation  requires  concepts  put  forth  by 
Brown  and  Lopez^  and  Darmofal.*®  Brown  and  Lopez^  established 
a  necessary  condition  between  the  production  of  negative  azimuthal 
vorticity  rj  and  the  extent  of  axial  flow  deceleration.  Consequently, 
the  authors  anticipated  that  the  minimum  azimuthal  vorticity  for 
the  three-dimensional  flow  would  be  greater  (less  negative)  than 
that  of  the  initial  two-dimensional  flow,  since  Q  is  greater  in  these 
cases.  This  is  indeed  the  case.  At  V  =  1.5,  before  the  Hopf  bi¬ 
furcation,  the  minimum  values  of  rj  computed  from  both  the  two- 
and  three-dimensional  models  are  virtually  the  same.  However,  just 
beyond  the  Hopf  point,  the  minimum  rj  is  significantly  greater  in 
the  three-dimensional  flow.  A  numerical  evaluation  of  asymmetric 
terms  in  the  azimuthal  vorticity  equation  is  then  ^rformed.^  This 
equation  relates  the  total  derivative  of  to  vorticity  stretching  and 
tilting  terms.*®  Through  this  evaluation,  it  is  found  that  the  effect  of 
three-dimensionality  on  radial  vorticity  is  the  principle  contributor 
to  the  increase  in  t).  A  region  consisting  of  a  positive  net  change  in 
radial  vorticity  exists  that  serves  to  attenuate  the  axial  deceleration 
process. 

Bowfield  visualization  is  performed  by  calculating  the  numerical 
equivalent  of  streaklines.  Bve  material  points  are  introduced  into  the 
flowfield  at  the  inflow  boundary.  The  white  point  lies  initially  on  the 
tube  centerline,  whereas  the  grey  scaled  points  lie  on  the  y  and  z  axes 
at  nondimensional  distances  of  ±0.1  from  the  centerline.  The  mate¬ 
rial  point  positions  are  computed  in  time  from  the  evolving  velocity 
field  using  a  first-order-accurate  Euler  time  integration.^  Snapshots 
of  the  time-asymptotic  streaklines  are  shown  in  Fig.  3;  the  tube 
geometry  is  omitted  for  clarity.  Figure  3  shows  the  development 
of  three-dimensional  flow  as  V  is  increased  past  the  Hopf  bifurca¬ 
tion  point.  At  V  =  2.3  (Fig.  3d),  the  material  points  deflect  off-axis 
in  a  well-defined  helical-type  structure,  consistent  with  spiral-type 
breakdown.  Further  discussion  of  the  flow  visualization  can  be  found 
in  Ref.  5. 
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Conclusions 

Time  integration  of  the  three-dimensional,  compressible  Navier- 
Stokes  equations  reveals  that  when  the  vortex  strength  is  increased 
past  a  critical  value,  the  time-asymptotic  flow  changes  from  ax¬ 
isymmetric  and  steady  to  asymmetric  and  time  periodic,  indicating 
a  supercritical  Hopf  bifurcation.  The  three-dimensional  flows  form 
a  solution  branch  that  bifurcates  from  the  path  of  two-dimensional 
solutions  at  the  bifurcation  point.  The  authors’  interpretation  of 
this  result  is  that  the  mechanism  for  the  existence  of  a  least  one 
family  of  three-dimensional  solutions,  which  possess  reversed  flow 
at  sufficiently  large  values  of  vortex  strength,  is  the  loss  of  sta¬ 
bility  of  the  axisymmetric  base  flowl  Minimum  values  of  axial 
velocity  Q  are  observed  to  increase  as  flow  asymmetries  develop 
just  beyond  the  Hopf  point.  Furthermore,  a  small  range  of  V  ex¬ 
ists  where  two-dimensional  solutions  exhibit  vortex  breakdown  but 
three-dimensional  solutions  do  not  This  attenuation  of  axial  decel¬ 
eration  is  found  by  the  authors  to  be  the  result  of  a  positive  net 
production  of  radial  vorticity  as  flow  asymmetries  develop. 
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The  three-dimensional,  compressible  Navier-Stokes  equations  in  primitive  variables  are  solved 
numerically  to  simulate  vortex  breakdown  in  a  constricted  tube.  Time  integration  is  performed  with 
an  implicit  Beam- Warming  algorithm  using  fourth-order  compact  operators  to  discretize  spatial 
derivatives.  Initial  conditions  are  obtained  by  solving  the  steady,  compressible,  and  axisymmetric 
form  of  the  Navier-Stokes  equations  with  Newton’s  method.  The  effects  of  three-dimensionality  on 
flows  that  are  initially  axisymmetric  and  stable  to  2-D  disturbances  are  examined.  Stability  of  the 
axisymmetric  base  flow  is  assessed  through  3-D  time  integration.  Axisymmetric  solutions  at  a  Mach 
number  of  0.3  and  a  Reynolds  number  of  1000  contain  a  region  of  nonuniqueness.  Within  this 
region,  3-D  time  integration  reveals  only  unique  solutions,  with  nonunique  axisymmetric  initial 
conditions  converging  to  a  unique  solution  that  is  steady  and  axisymmetric.  Past  the  primary  limit 
point,  which  approximately  identifies  the  appearance  of  critical  flow  (a  flow  that  can  support  an 
axisymmetric  standing  wave),  the  solutions  bifurcate  into  3-D  time-periodic  flows.  Thus  this 
numerical  study  shows  that  the  vortex  strength  associated  with  the  loss  of  stability  to  3-D 
disturbances  and  that  of  the  primary  limit  point  are  in  close  proximity.  Additional  numerical  and 
theoretical  studies  of  3-D  swirling  flows  are  needed  to  determine  the  impact  of  various  parameters 
on  dynamic  behavior.  For  example,  it  is  possible  that  a  different  flow  behavior,  leading  to  a  nearly 
axisymmetric  vortex  breakdown  state,  may  develop  with  other  inlet  profiles  and  tube  geometries. 
[81070-6631(97)02804-3] 


1.  INTRODUCTION 

Vortex  breakdown  is  a  hydrodynamic  feature  of  swirling 
flows  in  which  the  rotational  vortex  core  stagnates,  resulting 
in  a  dramatic  increase  in  the  core  size.  Breakdown  can  occur 
in  tornadoes,  in  swirling  flows  inside  tubes,  in  wing  trailing 
vortices,  and  in  the  vortical  flows  produced  over  delta  wings 
at  high  incidence. 

Many  numerical  studies  of  vortex  breakdown  rely  on  the 
assumption  of  a  steady  flow  which  is  symmetric  about  the 
vortex-core  axis.  This  simplification  allows  researchers  to 
study  the  nearly  axisymmetric  bubble  form  of  breakdown 
found  in  tube  experiments,^  while  avoiding  the  computa¬ 
tional  burden  of  solving  3-D  flowfields.  Some  of  the  more 
recent  studies  of  axisymmetric  breakdown  are  by  Buntine 
and  Saffman,^  Wang  and  Rusak,^’"^  Lopez, ^  Darmofal  and 
Murman,^  and  Beran  and  Culick,^  hereafter  denoted  as  BC. 

Three-dimensional  vortex  breakdown  computations  have 
also  been  reported.  These  computations  have  provided  valu¬ 
able  information  on  the  topological  structure  of  the  various 
3-D  breakdown  forms.  Most  notably,  calculations  have  been 
performed  for  flows  in  unconfined  domains^"* and  over 
delta  wings.  These  works  provide  descriptions  of  the  vari¬ 
ous  time-asymptotic  breakdown  structures  encountered,  such 
as  the  spiral  and  bubble  forms.  Spall  and  Gatski^^  compute 
other  flow  disturbances  in  addition  to  the  spiral  and  bubble 
breakdowns,  as  documented  in  the  tube  experiments  of  Faler 
and  Leibovich.^ 

An  interesting  feature  of  axisymmetric  swirling  flows  is 
the  appearance  of  nonunique  solutions  for  a  fixed  set  of  flow 
parameters.  Taasan^^  and  Leibovich  and  Kribus'^  computed 
solutions  to  the  Euler  equations  and  found  multiple  solutions 
stemming  from  bifurcations  of  columnar  flows.  BC  further 
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demonstrated  that  nonunique  flows  exist  as  solutions  of  the 
Navier-Stokes  equations.  These  results  also  indicate  a  con¬ 
nection  between  the  appearance  of  reversed  flow  and  the 
passage  of  a  limit  point.  The  stability  of  the  BC  solutions 
were  confirmed  by  Lopez^  and  Beran^"^  by  solving  the  time- 
dependent  form  of  the  2-D  governing  equations.  Lopez^  con¬ 
cluded  that  the  appearance  of  multiple  solutions  precludes 
the  possibility  of  predicting  breakdown  solely  from  knowl¬ 
edge  of  the  upstream  flow.  Rusak  and  Wang^^  recently  ap¬ 
plied  a  global  variational  analysis  to  incompressible,  inviscid 
flows.  They  proved  that  both  a  columnar  flow  and  a  flow 
with  a  large  stagnation  zone  can  coexist  within  a  specified 
range  of  vortex  strength.  These  two  solutions,  which  are 
stable  to  small  axisymmetric  disturbances,  are  connected  by 
an  unstable  branch  of  solutions. 

Nonuniqueness  in  3-D  swirling  flows  has  also  been  re¬ 
ported.  The  tube  experiments  of  Sarpkaya^^  show  hysteresis 
resulting  from  the  existence  of  both  the  bubble  and  spiral 
forms  of  breakdown  under  identical  flow  conditions.  Similar 
experiments  by  Faler  and  Leibovich^  show  spontaneous 
switching  between  the  spiral  and  bubble  forms  without  hys¬ 
teresis.  Nonuniqueness  in  delta-wing  flow  computations  is 
discussed  by  Visbal.^^ 

The  role  of  nonunique  axisymmetric  solutions  in  the  de¬ 
velopment  of  3-D  swirling  flows  is  the  subject  of  this  work. 
The  authors  have  previously  investigated  the  stability  of 
unique  axisymmetric  solutions  to  three-dimensional 
disturbances.^^  A  unique  feature  of  this  work  and  the  current 
work  is  that  solutions  to  the  steady,  axisymmetric  governing 
equations  are  used  as  initial  conditions  for  time-integration 
of  the  3-D  equations  of  motion.  Thus  transient  solutions 
which  deviate  away  from  the  initial  condition  do  so  as  a 
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Strength 

FIG.  1.  A  representative  solution  branch  of  axisymmetric  solutions;  stable 
branch  (solid  lines);  unstable  branch  (dashed  line). 


direct  consequence  of  three-dimensional  disturbances.  Time- 
asymptotic  3-D  solutions  are  compared,  under  identical  pa¬ 
rameter  settings,  to  the  2-D  initial  condition.  In  this  way,  the 
relevancy  of  2-D  solutions  is  revealed.  A  cylindrical  tube 
geometry  is  used,  since  it  should  naturally  give  rise  to  an 
axisymmetric  flow  for  a  sufficiently  low  vortex  strength. 
This  allows  for  a  clean  identification  of  the  emerging  vortex 
core  asymmetries  and  an  understanding  of  their  role  in  the 
development  of  breakdown. 

A  characterization  of  nonunique  solutions  for  viscous, 
steady  and  axisymmetric  flow  is  shown  in  Fig.  1,  plotted  in 
terms  of  the  (global)  minimum  axial  velocity  component, 
Q,  and  the  vortex  strength,  9^.  The  solid  lines  represent  so¬ 
lutions  that  are  stable  to  axisymmetric  disturbances,  while 
the  dashed  line  represents  unstable  (u)  solutions.  As  the  vor¬ 
tex  strength  is  increased,  a  fold  in  the  solution  space  occurs 
along  the  upper  stable  branch  (s*^),  resulting  in  a  primary 
limit  point  denoted  as  .  As  the  vortex  strength  is  slightly 
increased  from  Q  abruptly  changes  from  positive  to 
negative,  indicating  the  formation  of  bubble  breakdown.  For 
flows  with  breakdown  along  the  lower  stable  branch  (s'"), 
decreasing  the  vortex  strength  results  in  generally  larger  (less 
negative)  values  of  Q  until  the  secondary  limit  point,  5^^ ,  is 
encountered.  The  nonuniqueness  of  solution  paths  as  the  vor¬ 
tex  strength  is  increased  and  then  decreased  results  in  a  hys¬ 
teresis  loop.  Between  9^^  and  9^^ ,  nonunique  axisymmetric 
solutions  exist. 

The  primary  limit  point  in  an  axisymmetric  solution 
space  has  been  associated  with  the  appearance  of  critical 
flow — a  flow  state  which  supports  standing  axisymmetric 
waves,  A  supercritical  flow  therefore  refers  to  a  flow  which 
can  only  support  the  downstream  propagation  of  axisymmet¬ 
ric  waves,  while  subcritical  flows  can  support  either  up¬ 
stream  or  downstream  wave  propagation.  BC  found  that  for 
sufficiently  high  Reynolds  numbers,  a  parabolized  version  of 
the  axisymmetric  governing  equations,  known  as  the  quasi- 
cylindrical  (QC)  equations,  agree  well  with  solutions  of  the 
Navier-Stokes  equations  when  the  vortex  strength  is  below 
the  primary  limit  point.  However,  the  solutions  fail  to  con¬ 
verge  as  the  vortex  strength  is  increased  towards  the  primary 
limit  point.  The  failure  point  of  the  parabolic  QC  equations  is 
believed  by  Hall*^  to  be  the  approximate  point  at  which  criti¬ 
cal  flow  develops.  Thus,  the  works  of  BC  and  Hall'^  suggest 
that  the  flow  transitions  from  supercritical  flow  to  subcritical 
flow  at  a  vortex  strength  which  is  approximately  equal  to 
9^p .  More  recently,  a  noncolumnar  flow  criticality  analysis 
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has  been  reported  by  Wang  and  Rusak'^  for  2-D,  inviscid 
flow.  Their  results  indicate  that  the  critical  vortex  strength 
for  a  finite  length  tube  is  a  point  of  exchange  of  stability  for 
both  columnar  and  noncolumnar  solution  branches,  and  that 
this  critical  point  corresponds^  to  a  transcritical  bifurcation 
point  of  the  inviscid  solution  branches.  Furthermore,  an  ex¬ 
tension  of  their  work  to  flows  with  slight  viscosity^^  confirms 
the  existence  of  limit-point  behavior  found  previously  in  nu¬ 
merical  solutions  of  the  Navier-Stokes  equations.  The  pri¬ 
mary  limit  point  forms  at  a  slightly  lower  vortex  strength 
than  the  transcritical  bifurcation  point,  with  the  spacing  be¬ 
tween  the  two  determined  largely  by  the  magnitude  of  the 
viscosity.  The  primary  limit  point  is  also  expected to  be  a 
point  of  exchange  of  stability  to  axisymmetric  disturbances. 

The  computational  approach  of  the  current  study  is  as 
follows.  First,  a  numerical  algorithm  provides  the  axisym¬ 
metric  initial  condition  for  a  specified  vortex  strength,  9^’, 
given  fixed  values  of  the  freestream  Mach  number,  M,  and 
the  Reynolds  number.  Re.  The  axisymmetric  solution  is  then 
interpolated  onto  a  3-D  mesh  using  a  fourth-order-accurate 
cubic-spline  scheme.~^  Then,  time-integration  is  carried  out 
by  a  3-D  time-integration  algorithm.  Initial  conditions  are 
computed  using  the  pseudo-arclength  continuation  (PAC) 
model  of  BC  modified  for  compressible  flow.  The  PAC 
model  is  capable  of  computing  solution  paths  with  folds, 
providing  a  one-parameter  (vortex  strength)  family  of  axi¬ 
symmetric  solutions.  The  nominal  Mach  number,  M  =  0.3, 
was  chosen  to  provide  “near  incompressibility”  while 
avoiding  convergence  problems  associated  with  computing 
at  low  Mach  numbers.  The  3-D  time-integration  model,  re¬ 
ferred  to  as  the  Time-Accurate  Navier-Stokes  (TANS) 
model,  is  developed  specifically  for  this  work  and  is  de¬ 
scribed  in  the  next  section  and  in  the  Appendix. 

II.  NUMERICAL  MODEL 

The  TANS  model  incorporates  two  unique  features. 
First,  it  uses  a  multiblock  grid  structure  in  the  crossflow 
plane.  The  multiblock  structure  allows  for  a  nearly  Cartesian 
arrangement  of  nodes  in  the  vicinity  of  the  centerline,  while 
maintaining  near  orthogonality  at  the  circular  tube  wall.  This 
is  done  for  two  reasons:  (1)  to  allow  for  nearly  constant  grid 
resolution  near  the  centerline,  where  the  vortex  core  may 
migrate  off-center  during  spiral  breakdown;  and  (2)  to  avoid 
an  approximate  numerical  treatment  at  the  tube  centerline. 
Second,  the  TANS  model  incorporates  fourth-order  compact 
operators  into  an  approximate  factorization,  Beam-Warming 
solution  procedure.  The  compact  scheme  discretizes  explicit 
spatial  derivatives  to  fourth-order  accuracy.  Central- 
difference  discretizations  are  typically  second-order  accurate. 
The  net  reduction  in  CPU  time  using  the  compact  scheme 
over  central  differencing,  assuming  fewer  required  nodes  to 
achieve  similar  levels  of  accuracy,  is  about  42%. 

A.  Tube  geometry  and  grid 

A  Cartesian  (x,y,z)  and  a  cylindrical  (z,r,^)  coordinate 
system  are  referenced  such  that  positive  x  and  z  are  aligned 
with  the  tube  centerline  and  pointed  downstream.  The  y  and 
z  directions  lie  in  a  crossplane  normal  to  x  and  form  a  right- 
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of  clustering,  respectively.  Forx  =  0  andx>X2,  Ax,  =  Axo. 


FIG.  2.  Schematic  of  tube  geometry  in  an  arbitrary  (z,r)  plane. 


handed  system.  The  Cartesian  velocity  components  are  de¬ 
noted  by  y=(w,i;,w)^.  The  r  and  9  directions  are  oriented 
such  that  ^=0  corresponds  to  r=  +  y.  The  axial,  radial  and 
azimuthal  velocity  components  are  denoted  by  (vv',m,u),  re¬ 
spectively. 

The  physical  domain  consists  of  a  two-stage  cylindrical 
tube  of  circular  cross-section  and  varying  radius^  (Fig.  2). 
The  radius  of  the  inlet  station  is  denoted  as  /?o*  The  radius  of 
the  first  stage  is  given  by 

i?(x)  =  /?o+<^^o[cos(27r:c/xi)- 1]  (O^x^xi),  (1) 

where  Xi  is  the  length  of  the  first  stage.  The  parameter  a 
controls  the  amount  of  tube  contraction.  The  values  of  a  and 
Xi  are  fixed  for  this  work  at  0.05  and  6,18,  respectively.  The 
second  stage  of  the  tube  has  a  constant  tube  radius  of  Rq, 
The  domain  boundaries  are  denoted  by  si,  s2,  and  s3;  corre¬ 
sponding  to  the  inflow,  wall,  and  outflow  boundaries,  respec¬ 
tively  (Fig.  2).  A  generalized  mapping  transforms  the  physi¬ 
cal  coordinates  (jc, y,z)  to  the  computational  coordinates 
(|,77,f).  Node  indices  in  the  (^,77,^)  coordinates  are  de¬ 
noted  by  respectively.  The  number  of  nodes  in  the 

(^,77,^)  coordinate  directions  are  (nx,ny,nz),  respectively. 

A  crossplane  of  the  multiblock  grid  is  shown  in  Fig.  3(a) 
for  «y  =  «z  =  41  .  The  grid  consists  of  an  inner  block  sur¬ 
rounded  by  four  outer  blocks  [Fig.  3(b)].  The  outer  blocks 
are  physically  connected  to  each  other  but  contain  two  edges 
each  that  are  considered  as  branch  cuts  in  the  computational 
domain.  These  cuts  are  labeled  1-4  in  Fig.  3(b).  The  2-D 
crossplane  grid  is  generated  with  the  GRIDGEN  (Ref.  21)  soft¬ 
ware  package.  The  3-D  grids  are  generated  by  scaling  the 
crossplane  grid  to  the  local  tube  radius.  Clustering  is  allowed 
in  the  x  direction.  The  axial  node  locations  are  computed 
from  1 4-  Ax,- ,  where 


FIG.  3.  Crossplane  grid:  (a)  physical,  (b)  computational. 


B.  Governing  equations 

The  nondimensional  vector  form  of  the  Navier-Stokes 
equations  is  written  below.  The  fluid  density  and  pressure  are 
denoted  by  p  and  p,  respectively,  to  define  the  solution  vec¬ 
tor  f/=(p,w,i;,w,p)^.  For  algorithm  efficiency,  the  equa¬ 
tions  are  written  in  nonconservative  form,  limiting  the 
scheme  to  subsonic  flow.  The  reference  length  is  the  inlet 
radius  of  the  vortex  core.  The  reference  velocity  is  the  (uni¬ 
form)  inflow  axial  velocity  and  the  reference  pressure  is 
twice  the  reference  dynamic  pressure 
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where  is  the  Kronecker  delta  function  for  /,m=  1,2,3  and 
D  =  [0,V  •  T,p{(  7—  !)<!>  — V  •q}y.  The  nondimensional 
shear  stress  tensor,  r,  is  defined  assuming  Stokes’  hypoth¬ 
esis,  $  and  q  represent  the  nondimensional  viscous  dissipa¬ 
tion  and  heat  flux  vector,  respectively.  The  auxiliary  equa¬ 
tions  necessary  to  close  the  system  of  equations  come  from 
the  assumption  of  a  perfect  gas  and  Sutherland’s  formula. 

Axisymmetric  conditions  are  specified  at  the  inflow 
boundary  plane,  si.  The  inlet  axial  velocity,  m  =  is  chosen 
to  be  uniform.  Boundary  conditions  on  v  and  w  are  obtained 
by  specifying  appropriate  profiles  for  v  and  u.  The  following 
swirl  velocity  profile  is  characteristic  of  a  Burger  type  vortex 
and  is  appropriate  for  modeling  the  profiles  obtained  from  a 
swirl- vane  apparatus:^ 

where  W'  is  the  vortex  strength  along  si  and  F  is  the  circu¬ 
lation  (divided  by  27r).  The  radial  velocity,  w,  is  assumed  to 
vanish  to  reflect  a  columnar  flow  state.  A  nonuniform  inflow 
density  profile  based  on  columnar  flow  is  prescribed  as  a 
Dirichlet  condition.  The  profile  is  obtained  by  solving  the 
axisymmetric  Navier-Stokes  equations  with  the  PAC  algo¬ 
rithm  in  a  straight  tube  of  short  axial  extent  (L  —  0.01).  In¬ 
flow  and  outflow  conditions  are  enforced  that  dictate  colum¬ 
nar  flow.  These  conditions  consist  of  fixed  velocity 
components  (given  above)  with  py.=p^=().  The  resulting  co¬ 
lumnar  solution  for  density  is  then  used  as  a  Dirichlet  bound¬ 
ary  condition.  A  boundary  condition  on  pressure  is  obtained 
by  solving  the  steady  form  of  the  axial  momentum  equation 
along  the  inflow  plane.  The  constriction  in  the  first  stage  of 
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the  tube  provides  a  favorable  pressure  gradient  in  the  con¬ 
verging  section,  while  an  adverse  gradient  forms  aft  of  the 
throat.  Thus  the  constriction  keeps  the  breakdown  from  oc¬ 
curring  near  si,  where  columnar  conditions  are  enforced.  A 
second  inlet  condition  of  zero  azimuthal  vorticity  is  briefly 
considered,  implying  a  vanishing  axial  derivative  of  radial 
velocity  for  a  uniform  axial  velocity  profile.  This  condition 
relaxes  the  columnar  flow  assumption,  allowing  for  a  com¬ 
parison  of  results  between  two  inflow  conditions.  The  colum¬ 
nar  conditions  are  assumed  unless  otherwise  noted.  Surface 
s2  is  assumed  to  be  a  slip-wall,  yielding  an  impermeability 
boundary  condition  which  relates  the  axial  and  radial  veloci¬ 
ties.  Two  other  conditions,  based  on  the  compressible  Ber¬ 
noulli  equation  and  the  assumption  of  constant  wall  circula¬ 
tion,  lead  to  three  Dirichlet  conditions  for  the  velocity 
components  at  the  wall.  A  fixed  wall  temperature  is  also 
assumed.  The  wall  pressure  is  determined  by  solving  the 
steady  form  of  the  y  and  z  momentum  equations.  The  y 
momentum  equation  is  solved  in  the  77  sweep  and  the  z  mo¬ 
mentum  equation  is  solved  in  the  ^  sweep.  Solution  sweeps 
are  discussed  in  the  Appendix.  Outflow  conditions  at  s3'  are 
chosen  to  reflect  an  assumed  columnar  flow  state.  Details  of 
the  boundary  condition  formulation  can  be  found  in  Ref.  20. 

111.  RESULTS 

The  following  results  show  that  steady  solutions  occur 
when  the  vortex  strength  is  prescribed  to  be  less  than  the 
primary  limit  point.  In  particular,  the  specification  of  three 
nonunique,  axisymmetric  initial  conditions  (at  a  vortex 
strength  between  the  primary  and  secondary  limit  points) 
lead  to  three  apparently  identical  time-asymptotic  solutions 
which  are  steady  and  axisymmetric.  When  the  vortex 
strength  is  prescribed  to  be  greater  than  the  primary  limit 
point,  a  time-periodic,  three-dimensional  solution  develops. 
The  change  in  solution  behavior  is  attributed  to  a  Hopf  bi¬ 
furcation,  which  is  found  to  lie  in  very  close  proximity  to  the 
primary  limit  point. 

A  Reynolds  number  of  1000  is  considered  following  the 
results  of  BC,  who  show  that  nonunique  2-D  solutions  exist 
for  incompressible  flow  when  Re>360.  Consistent  results 
are  found  for  M  =  0.3.  Only  unique  axisymmetric  solutions 
were  evident  in  previous  results  at  Re=250,^^  while  a  region 
of  nonunique  solutions  are  found  at  Re =500  (Ref.  20)  and 
Re=  1000.  The  vortex  strength  difference  between  primary 
and  secondary  limit  points  is  0.0019  and  0.015  for 
Re  =  500  and  Re=  1000,  respectively.  BC  report  limit  point 
separations  for  M  =  0  of  about  0.012  and  0.04  for  Re  =500 
and  Re  =1000,  respectively,  which  are  greater  than  those 
found  here  for  Af  =  0.3.  Thus  compressibility  effects  appear 
to  diminish  the  region  of  nonuniqueness  for  a  given  Rey¬ 
nolds  number. 

A,  Grid  requirements 

Grid  requirements  for  2-D  swirling  flows  have  been  dis¬ 
cussed  in  previous  work  by  BC,  Lopez, ^  and  Darmofal  and 
Murman.^  Each  of  these  works,  along  with  the  current  study, 
consider  the  constricted  tube  geometry  defined  by  BC.  More 
specifically,  the  outer  tube  radius  is  governed  by  Eq.  (1), 
a=0.05  and  Rq-I.  Darmofal  and  Murman^  also  consider  a 


second  tube  constriction  near  the  outflow  boundary.  BC 
demonstrate  solution  insensitivity  by  prescribing  27X301 
nodes  for  Ro  =  2  and  L  =  30.  Some  sensitivity  between  tube 
lengths  L  =  30  and  L  =  60  are  noted  at  Re  =2000,  which  is 
higher  than  the  Re=  1000  flows  considered  here.  Similarly, 
Lopez^  used  31X301  nodes  for  Rq-2  and  L  =  30.  Further 
resolution  of  this  grid  to  61X601  removed  the  appearance  of 
spatially  developing  “wiggles”  in  contours  of  azimuthal 
vorticity.  However,  Lopez^  points  out  that  the  differences 
between  the  31X301  solutions  and  the  61X601  solutions  are 
not  significant  for  Re  =^1000.  In  particular,  plots  of  the 
streamfunction  and  circulation  for  the  two  grids  were  de¬ 
scribed  as  being  virtually  indistinguishable.  BC  also  found 
that  reducing  the  number  of  axial  nodes  from  301  to  151 
axial  nodes  had  very  little  effect  on  the  solution  paths,  except 
that  solution  convergence  suffered  near  the  secondary  limit 
point.  Darmofal  and  Murman^  used  31X151  nodes  for 
Rq  =  2  and  L  =  30  to  compute  solutions  at  Re=1000. 

Grid  requirements  are  obviously  greater  for  the  current 
3-D  flows.  However,  using  the  fourth-order  compact  opera¬ 
tors  provides  some  relief  by  allowing  similar  accuracy  levels 
with  fewer  nodes.  In  this  study  a  fine  grid  is  constructed  with 
61X61X172  nodes  for  Rq  =  2  and  L=30.  This  grid  is 
equivalent  in  crossplane  resolution  to  the  above  grids.  Grid 
clustering  in  the  axial  direction  is  employed  [Aa:o  =  0.206, 
(3  =  0.25y  and  jc2=10;  Eq.  (2)]  near  the  aft  portion  of  the 
throat  such  that  the  spacing  at  x  =  5  is  one-half  the  spacing 
achieved  with  equally  spaced  nodes.  Thus  the  minimum  grid 
spacing  in  the  breakdown  region  is  comparable  to  the  spac¬ 
ing  achieved  with  301  equally  spaced  nodes.  Using  the 
fourth-order  compact  operators  further  enhances  accuracy. 
Results  in  Ref.  20  show  that  this  level  of  axial  resolution  is 
adequate  for  Re=^  1000.  A  time  step  of  Af  =  0.025  is  used, 
based  on  time-accuracy  requirements  identified  in  Ref.  20. 

Varying  the  tube  length  from  L  =  20  to  L  =  30  results  in 
large  differences  in  Q  when  the  vortex,  strength  is  greater 
than  the  primary  limit  point  and  when  Re=  1000.  However,  a 
tube  length  of  L  =  20  was  found  to  be  adequate  for 
Re=250,^^  Re=500,^®  and  for  sufficiently  low  values  of 
vortex  strength  at  Re=  1000.  The  observed  sensitivity  to  tube 
length  prompted  all  runs  computed  for  Re=  1000  to  be  run 
with  L  =  30.  A  further  extension  of  the  tube  to  L  =  40  shows 
that  only  small  deviations  in  Q  occur  from  the  L  =  30 
solution. 

Values  of  tube  radius  other  than  Rq  =  2  are  not  consid¬ 
ered  in  this  study.  BC  found  that  paths  of  axisymmetric  so¬ 
lutions  near  the  primary  limit  point  are  insensitive  to  tube 
radius  between  Ro  =  '^  and  Ro  =  3  for  Reynolds  numbers 
ranging  from  250  to  2000.  However,  some  sensitivity  in  the 
bifurcation  point  locations  of  inviscid,  axisymmetric  solution 
branches  to  the  vortex  core/tube  radius  ratio  have  been 
reported.^^  At  vortex  strengths  much  higher  than  the  primary 
limit  point  (^ as  high  as  2.3),  3-D  flows  were  computed^®  at 
Re=250,  without  the  anticipated  formation  of  the  bubble 
form  of  breakdown.^  Bubble  breakdowns  were  instead  only 
observed  as  transient  solutions,  and  ultimately  transformed 
to  one  of  the  asymmetric  (spiral)  modes.  Evidence  of  the 
transient  bubble  form  suggests  that  the  level  of  grid  resolu¬ 
tion  is  adequate  to  resolve  the  bubble  breakdown.  Instead, 
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FIG.  4.  Q  versus  time:  (dashed),  (solid).  9''^p  —  1.47435. 


the  authors  speculate  that  the  relatively  small  value  of  the 
tube  radius,  Rq  —  2,  is  insufficient  to  produce  a  3-D  time- 
asymptotic  bubble  form  of  breakdown  at  high  vortex 
strengths.  The  basis  for  this  belief  is  the  observation  that 
flowfield  perturbations  near  the  tube  wall  are  larger  for  3-D 
solutions  than  for  2-D  solutions.  Thus,  further  computations 
with  the  model  should  consider  larger  values  of  tube  radius 
to  determine  if  the  wall  is  preventing  the  occurrence  of  the 
3-D  bubble  breakdown  at  large  values  of  It  is  of  interest 
to  note  that  a  small,  but  steady,  axisymmetric  bubble  was 
computed  with  the  TANS  model  for  a  much  lower  Reynolds 
number  of  50. 

B.  Temporal  behavior 

Six  axisymmetric  initial  conditions  are  time-integrated 
with  the  TANS  model  for  Re=  1000.  Three  of  these  initial 
conditions  lie  before  the  primary  limit  point  at 
^p==  1.47435.  These  nonunique  solutions  lie  on  the  upper 
stable  branch  1.47'^'^),  the  unstable  branch 

1.47“),  and  the  lower  stable  branch  1.47"-').  The 
remaining  solutions  are  obtained  at  ^=1.48,  9^=1.5,  and 
1.6,  which  are  all  greater  than  . 

The  temporal  behavior  of  Q  for  9^=1.47  and  1.48  is 
shown  in  Fig.  4.  It  is  evident  that  the  flow  temporally  evolves 
from  a  negative  to  a  positive  value  of  Q  for  9^=IAT~, 
9^ =  1.47“  and  9^ =  1.48,  indicating  the  elimination  of  a  re¬ 
versed  flow  region.  Furthermore,  the  three  runs  computed 
before  the  primary  limit  point  (dashed  lines  in  Fig.  4)  ail 
migrate  to  apparently  identical  solutions  near  the  initial  con¬ 
dition  at  9^=  1.47^'‘'.  This  behavior  signifies  the  loss  of  2-D 
solution  nonuniqueness  in  the  presence  of  three-dimensional 
disturbances.  When  9^  is  increased  slightly  to  1,48,  the  tem¬ 
poral  behavior  of  Q  departs  from  the  behavior  at 
9^=IA1^~  near  /  =  70,  The  time-asymptotic  value  of  Q 
(0.195  at  r=1000)  for  9^— IAS  is  lower  than  the  time- 
asymptotic  values  computed  before  the  primary  limit  point  at 
9^=  1.47,  but  is  higher  than  the  initial  2-D  value.  Thus,  the 
initially  steady,  axisymmetric  bubble  breakdown  solution  is 
replaced  by  a  flow  without  breakdown  when  three- 
dimensionality  is  allowed. 

The  dynamics  of  the  9^—lAT~  case  is  examined  to 
better  understand  the  evolution  from  axisymmetric  bubble 
breakdown  to  nearly  columnar  flow.  The  numerical  equiva¬ 
lent  of  streaklines  are  computed  for  this  run  to  provide  a 
visual  representation  of  the  flowfield.  Five  “material  points’" 
are  introduced  into  the  flowfield  at  the  inflow  boundary.  One 
point  (black)  lies  initially  on  the  tube  centerline,  the  others 
(gray  scaled)  lie  on  the  y  and  z  axes  at  a  fixed  offset  from  the 
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FIG.  5.  Streaklines  showing  transition  from  axisymmetric  bubble  break¬ 
down  to  nearly  columnar  flow:  1.47^“,  Re=  1000. 


centerline.  Five  new  material  points  are  periodically  intro¬ 
duced  into  the  flowfield  where  the  first  five  points  were  ini¬ 
tialized,  thereby  simulating  a  numerical  equivalent  of  die- 
injection  in  experiments.  Material  point  positions  are 
computed  in  time  from  the  evolving  velocity  field  using  a 
first-order-accurate  Euler  time  integration.^®  Initial  streak- 
lines  (r=0)  are  computed  from  the  initially  steady  flowfield. 
When  the  material  points  are  initially  offset  from  the  center- 
line  a  nondimensional  distance  of  0.1,  the  resulting  streak¬ 
lines  convect  over  the  bubble  breakdown  region  and  down¬ 
stream.  However,  if  the  initial  offset  is  0.01,  material  points 
pass  over  the  bubble  and  then  enter  the  bubble  region  from 
the  rear.  Both  offsets  are  considered  to  allow  material  point 
visualizations  inside  and  outside  the  bubble. 

Examination  of  the  streaklines  at  discrete  points  in  time 
reveal  that  the  initially  axisymmetric  bubble  breakdown 
erupts  temporarily  into  a  3-D  spiral  flow,  followed  by  a  de¬ 
cay  of  three-dimensionality  to  the  nearly  columnar  solution 
along  the  s*^  branch.  Furthermore,  the  instability  associated 
with  the  appearance  of  3-D  flow  appears  greatest  just  aft  of 
the  breakdown  bubble. 

The  emergence  and  decay  of  3-D  flow  is  illustrated  in 
Fig.  5,  At  r=0,  the  bubble  is  clearly  evident  by  streaklines 
which  pass  around  the  bubble  region  and  convect  down¬ 
stream.  The  bubble  structure  remains  intact  through  r  =  25. 
However,  as  discussed  later,  it  is  found  that  asymmetries 
have  developed  both  within  and  aft  of  the  bubble  region  by 
t-25  and  have  significantly  altered  the  flow  inside  the 
bubble.  Thus  the  material  points  shown  in  Fig.  5  do  not 
provide  a  total  picture  of  the  emerging  3-D  flow.  By  r  =  40, 
the  aft  portion  of  the  bubble  is  replaced  with  a  loosely  orga¬ 
nized  spiral  arrangement  of  material  points.  The  material 
points  upstream  and  slightly  downstream  of  the  nose  of  the 
bubble  remain  nearly  undisturbed.  By  t  =  50,  a  spiral  struc¬ 
ture  replaces  the  bubble  entirely,  and  then  slowly  transitions 
over  time  to  nearly  columnar  flow.  By  r=60,  the  original 
location  of  the  nose  of  the  bubble  is  replaced  with  nearly 
columnar  flow.  The  axial  extent  of  the  nearly  columnar  flow 
increases  downstream  with  time  as  evidenced  by  the  streak¬ 
lines  at  100.  The  flow  becomes  steady  and  nearly  colum¬ 
nar  everywhere  by  r  =  400. 
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FIG.  6.  Centerline  values  of  the  y  velocity  component,  u,  versus  x  for 
various  times,  t\  1.47^",  Re=  1000. 

The  emergence  of  3-D  flow  aft  of  the  bubble  is  illus¬ 
trated  in  Fig.  6,  where  the  centerline  value  of  the  y  velocity 
component,  u,  is  plotted  versus  .r  at  various  times.  Note  that 
nonzero  centerline  values  of  v  imply  3-D  flow.  The  location 
of  the  breakdown  bubble  is  marked  in  the  figure  by  two 
vertical  dashed  lines.  It  is  found  that  below  r  =  5,  small 
asymmetries,  exist  slightly  upstream  and 

downstream  of  the  bubble  region.  As  t  is  increased,  the 
asymmetries  upstream  of  the  bubble  grow  to  only 
Within  and  especially  aft  of  the  bubble  region, 
however,  asymmetric  flow  develops  quickly.  By  f  =  25, 
asymmetric  flow  extends  from  the  front  portion  of  the  bubble 
region  to  as  far  downstream  as  ;c  =  20.  This  asymmetric  flow 
develops  before  any  obvious  signs  of  their  existence  appears 
in  the  streaklines  of  Fig.  5. 

The  streaklines  shown  in  Fig.  7  convect  upstream  into 
the  interior  of  the  bubble,  and  help  visualize  the  impact  of 
the  emerging  asymmetric  flow.  Between  r  =  0  and  t=5  the 
bubble  structure  remains  basically  unchanged  from  the  initial 
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FIG.  7.  Streaklines  showing  flow  reversal  into  the  breakdown  bubble  and 
subsequent  affect  of  emerging  flow  asymmetries:  7'*=  1.47^",  Re=  1000. 
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condition,  with  all  material  points  entering  the  bubble  from 
the  rear  and  spiraling  upstream.  By  r=  15,  two  changes  to  the 
aft  portion  of  the  bubble  are  apparent.  First,  the  spiraling 
material  points  within  the  bubble  between  approximately 
x  =  S  and  x  =  9  are  not  as  tightly  wrapped  as  before.  This  is 
found  to  be  associated  with  the  elimination  of  reversed  flow 
near  the  rear  of  the  bubble.  Second,  the  rear  of  the  bubble 
appears  to  have  shifted  slightly  downstream.  By  r  =  20,  re¬ 
versed  flow  is  lost  in  the  downstream  half  of  the  bubble 
interior,  and  the  ejection  of  material  points  out  of  the  bubble 
interior  is  evident.  Furthermore,  three-dimensional  flow  is 
now  clearly  evident  from  the  material  points.  By  f  =  25,  the 
region  of  reversed  flow  along  the  centerline  is  contained  to 
the  region  between  the  nose  of  the  bubble  and  x^l. 

The  dynamic  behavior  illustrated  above  in  Figs.  5-7 
demonstrates  the  instability  of  the  initial  axisymmetric  flow 
along  the  s~  branch  to  3-D  disturbances.  Furthermore,  the 
instability  is  strongest  just  aft  of  the  breakdown  bubble  (Fig. 
6). 

The  explanation  for  the  subsequent  decay  of  the  3-D 
disturbances,  however,  is  not  clear.  A  possible  explanation  is 
as  follows.  A  large  increase  in  Q  is  observed  between 
r^40  and  f*=^50  (Fig.  4),  and  is  found  to  coincide  with  the 
largest  growth  of  flow  asymmetries.  The  increase  in  Q  may 
be  a  direct  result  of  the  emerging  3-D  flow,  as  was  found  in 
previous  work  by  the  authors. As  a  result  of  increasing 
Q,  the  flow  may  reach  a  point  in  time  where  the  underlying 
axisymmetric  base  flow  is  no  longer  unstable  to  3-D  distur¬ 
bances.  [It  is  of  interest  to  note  that  flow  asymmetries  begin 
to  diminish  at  a  point  in  time  (r^47)  when  Q  switches  sign 
from  negative  to  positive.]  This  would  then  cause  the  subse¬ 
quent  decay  of  the  spiral  formation  in  Fig.  5  to  axisymmetric 
flow.  However,  the  strongest  attracting  axisymmetric  solu¬ 
tion  is  now  the  nearly  columnar  solution  along  the  s*^ 
branch.  In  summary,  solutions  along  s'"  are  unstable  to  3-D 
disturbances.  The  emerging  asymmetries  are  associated  with 
an  increase  in  Q  (i.e.,  away  from  the  s“  branch)  until  a 
sufficiently  high  value  of  Q  is  reached.  At  this  value  of  <2  >  it 
is  postulated  that  the  base  flow  can  no  longer  support  asym¬ 
metric  disturbances.  The  solution  would  then  be  attracted  to 
the  nearly  columnar,  s"^,  branch. 

A  similar  scenario  may  exist  for  a  flow  computed  along 
the  s“  branch  past  the  primary  limit  point.  In  this  case,  how¬ 
ever,  there  is  no  attracting  s’**  solution  available,  providing 
the  scenario  for  a  new  branch  of  solutions  to  form. 

C.  Evidence  for  a  Hopf  bifurcation  near  the  primary 
limit  point 

The  time-asymptotic  behavior  of  the  flows  computed 
near  the  primary  limit  point  is  discussed  here  and  compared 
to  the  initial  flow  state.  The  time-asymptotic  behavior  of  the 
flows  computed  before  the  primary  limit  point  is 

steady  and  axisymmetric.  Flow  steadiness  is  confirmed  by 
monitoring  the  global  maximum  value  of  the  solution  incre¬ 
ment,  A”t/,  as  time  increases.  All  of  the  runs  computed  be¬ 
fore  the  primary  limit  point  converge  to  maximum  values  of 
A"f/  near  1  X  10“^  at  r  =  400,  which  is  considered  an  accept¬ 
able  level  of  convergence.  The  time-asymptotic  solutions  are 
found  to  be  axisymmetric  by  computing  the  value  of  //, 
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FIG.  8.  Periodic  behavior  of  flow  computed  at  ^^=1.5  and  Re=  1000:  (a) 
Uj  (solid)  and  w,  (dashed)  versus  time,  (b)  phase  plot  of  for  r  =  370  to 
r  =  400. 


which  characterizes  the  degree  of  flow  asymmetry  and  is 
defined  as  the  global  maximum  absolute  value  of  dvIdO, 
H  is  found  to  be  very  close  to  zero  when  the  vortex  strength 
is  below  For  example,  for  the  time- 

asymptotic  value  of  //  is  8.4X  Furthermore,  it  is  found 
that  89%  of  this  value^®  is  due  to  interpolation  errors  in  com¬ 
puting  H.  The  remaining  11%  is  believed  to  be  attributable 
to  the  asymmetric  arrangement  of  the  grid  nodes,  and  not  to 
the  physical  presence  of  asymmetric  flow. 

The  time-asymptotic  behavior  past  the  primary  limit 
point  is  three-dimensional  and  time-periodic.  Values  of  H  for 
time-asymptotic  flows  computed  past  the  primary  limit  point 
are  much  greater  than  the  near-zero  value  of  H  computed  at 
9^=1.47""'.  For  example,  //=  1.293X  10” ^  for  the  time- 
asymptotic  flow  at  1.48.  The  time-periodic  nature  of  the 
flow  at  is  observed  by  plotting  time  histories  of  flow  quanti¬ 
ties  at  a  fixed  location  in  space.  In  particular,  Vs  and  are 
defined  as  the  crossplane  velocity  components  at  a  fixed  cen¬ 
terline  location  of  x  =  7.  Plots  of  and  are  shown  in  Fig. 
8(a),  with  the  corresponding  phase  plot  of  shown  in  Fig. 
8(b)  (<5r=0.8).  The  periodic  nature  of  is  evident  in  the 
figure.  The  phase  plot  of  is  not  shown,  owing  to  the 
similarity  between  and  .  The  axial  velocity  component 
Wj  is  found  to  be  steady. 

Contours  of  the  axial  velocity  component,  w,  are  shown 
in  Fig.  9  to  illustrate  the  change  in  flow  behavior.  Contours 
of  u  are  shown  in  Fig.  9(a)  for  9^=  lAT'^.  One-half  of  the 
x-y  plane  is  omitted  due  to  flow  symmetry  about  the  a:  axis. 
The  initial,  steady  axisymmetric  flow  is  shown  with  dashed 
contour  lines,  while  the  time-asymptotic,  steady  solution  of 
the  TANS  model  is  shown  with  solid  contour  lines.  The 
agreement  between  the  two  solutions  is  excellent,  demon¬ 
strating  that  the  initial  flow  is  stable  to  the  small  3-D  distur¬ 
bances  introduced  in  the  numerical  solution  procedure.  Fig. 
9(b)  shows  the  initial  axisymmetric  solution  for  ^=1.5, 
which  is  past  the  primary  limit  point.  The  strong  bubble 
breakdown  region  is  clearly  evident  from  the  negative  region 
of  u  (dashed  lines)  between  x  =  5  and  ;c=  10.  Time  integra¬ 
tion  of  this  flowfield  to  r  =  400  results  in  the  3-D  flowfield 
depicted  in  Fig.  9(c).  The  bubble  breakdown  structure  van¬ 
ishes  when  3-D  flow  is  allowed,  resulting  in  a  solution  with 
no  reversed  flow.  The  three-dimensional  nature  of  the  flow  is 
found  to  be  contained  to  the  region  aft  of  the  tube  constric- 
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FIG.  9.  Contours  of  the  axial  velocity  component,  w,  in  the  x-y  plane 
(z~0)  for  Re=  1000  flows  computed  about  1.47435  (numbers  in  pa¬ 
rentheses  indicate  the  range  of  15  evenly-spaced  contour  values):  (a) 
9''~\AT'^  (0.44,1.4);  initial  axisymmetric  solution  (dashed),  TANS  time- 
asymptotic  solution  (solid),  (b)  initial  axisymmetric  solution  at  1.5 
(-0.12,1,4);  w>0  (solid),  w<0  (dashed),  (c)  3-D  TANS  time-asymptotic 
solution  for  1.5  at  r=400  (0.1, 1.4). 


tion,  with  steady  and  axisymmetric  flow  maintained  up¬ 
stream. 

Snapshots  of  the  time-asymptotic  streaklines  for 

1.47^‘’'  and  1.5  are  shown  in  Fig.  10;  the  tube  ge¬ 
ometry  is  omitted  for  clarity.  The  flows  have  a  base  rotation 
in  the  clockwise  direction  when  viewed  in  the  positive  x 
direction.  At  ^=1.47'"'^  [Fig.  10(a)],  before  the  Hopf  point, 
the  vortex  core  swells  at  an  axial  location  of  about  x= 6.  The 
swelling  of  the  core  occurs  symmetrically,  and  illustrates  the 
effect  of  axial  flow  deceleration.  At  9^'=1.5  [Fig.  10(b)], 
beyond  the  Hopf  point,  a  spiral  asymmetric  disturbance  oc¬ 
curs  upstream  of  the  initial  swelling  at  x=6  for 
^=1.47^'^.  The  upstream  movement  of  flow  disturbances 
with  increasing  evident  between  Figs.  10(a)  and  10(b),  is 
in  agreement  with  experimental  evidence.^  The  axial  decel¬ 
eration  along  with  a  time-periodic  rotation  of  the  disturbance 
produces  distorted  rings  of  material  points,  which  subse¬ 
quently  spiral  and  convect  downstream.  The  direction  of  ro¬ 
tation  is  clockwise  (looking  downstream);  consistent  with  the 


FIG.  10.  Particle  traces  showing  development  of  three-dimensionality  be¬ 
yond  the  primary  limit  point  at  1.47435:  (a)  1.47''‘^,  (b)  1.5. 

J.  C.  Tromp  and  P.  S.  Beran 


1.46  1.48  1.50  1.3  1.4  1.5  1.6  1.7 

Vortex  Strength  Vortex  Strength 

FIG.  ll.  Bifurcation  diagram:  (a)  Initial  and  time-asymptotic  values  of  Q 
showing  emergence  of  3-D  solution  branch  (diamonds)  from  branch  of  axi- 
symmetric  solutions  (squares),  (b)  magnitude  of  time-periodic  disturbances, 
n,  (squares)  and  magnitude  of  flow  asymmetries,  H  (circles). 

base  vortical  flow.  Solutions  at  higher  vortex  strengths  reveal 
more  coherent  spiral  disturbances,  similar  to  previous  results 
at  Re=250.'’ 

The  loss  of  stability  of  steady  flow  to  time-periodic  flow 
as  vortex  strength  is  increased  is  evidence  for  a  Hopf  bifur¬ 
cation.  The  nature  of  the  bifurcation  is  summarized  in  Fig. 
11.  In  Fig.  11(a),  the  initial  and  time-asymptotic  values  of 
Q  are  depicted.  The  initial  conditions  are  denoted  with 
square  symbols  while  time-asymptotic  values  are  time- 
averaged  and  denoted  with  diamond  symbols.  A  distinct 
change  in  solution  character  occurs  when  9^^  is  increased 
slightly  from  1.47  to  1.48.  All  three  time-asymptotic  solu¬ 
tions  at  9"  =  1.47  are  steady  and  axisymmetric,  while  at 
^^’=1.48  the  time-asymptotic  flowfield  is  time-periodic  and 
three-dimensional.  In  between  the  vortex  strengths  of  1.47 
and  1.48  lies  the  primary  limit  point  at  1.47435.  Thus 
the  maximum  difference  in  9'"  between  the  location  of  the 
Hopf  bifurcation  and  the  primary  limit  point  is  5.65 X  10“^. 
The  diamond  symbols  denote  the  emerging  path  of  3-D  so¬ 
lutions  which  bifurcate  from  the  branch  of  steady,  axisym¬ 
metric  solutions.  The  authors’  interpretation  of  this  result  is 
that  the  mechanism  for  the  existence  of  at  least  one  family  of 
3-D  solutions,  which  possess  reversed  flow  at  sufficiently 
large  values  of  vortex  strength,  is  the  loss  of  stability  of  the 
axisymmetric  base  flows.  Reversed  flow  is  not  observed 
along  the  3-D  solution  branch  until  the  vortex  strength  is 
increased  from  5^''=  1.5  to  9^"'=  1.6.  After  the  primary  limit 
point,  ^=1.48,  the  time-asymptotic  value  of  Q  (diamond 
symbol)  is  significantly  larger  than  the  initial  value  of  Q 
(square  symbol)  for  the  same  value  of  9'\  This  observation  is 
consistent  with  previous  results^^  where  the  increase  in  Q  is 
correlated  to  the  positive  net  production  of  both  radial  and 
azimuthal  vorticity  components  as  flow  asymmetries  emerge. 

The  nature  of  the  Hopf  bifurcation  is  further  illustrated 
by  plotting  the  parameters  H  and  D  in  Fig.  11(b).  Cl  is 
defined  as  the  maximum  and  minimum  centerline  values  of 
the  sum  u  +  vw  found  over  a  sufficiently  long  time  interval  in 
the  time-asymptotic  flow.  Thus  fl  is  zero  for  axisymmetric 
flow  and  nonzero  for  3-D  flow.  The  fact  that  D  steadily 
increases  from  zero  near  9^  is  evidence  for  the  supercritical 
type  of  Hopf  bifurcation.^^  Furthermore,  it  is  evident  that 
H  sharply  departs  from  zero  for  flows  computed  past  the 
primary  limit  point.  Values  of  //,  D  and  Q  for  Re=  1000  are 
summarized  in  Table  I. 

The  solid  line  in  Fig.  11(b)  is  a  quadratic  polynomial  fit 
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TABLE  I.  Summary  of  time-asymptotic  parameter  values  for  Re=  1000. 


■r 

Q  (2-D)/(3-D) 

n 

H 

\Ar* 

-h0.369/+0.371 

-0.0001/0.0001 

0.0084 

1.47“ 

+0.056/+0.364 

-0.0001/0.0002 

0.0084 

1.47^" 

-0.098/+0.369 

-0.0003/0.0004 

0.0084 

1.48 

-0.110/+0.195 

-0.1550/0.1550 

0.1293 

1.50 

-0.128/+0.124 

-0.2513/0.2513 

0.2491 

1.60 

-0.159/-0.032 

-0.5080/0.5080 

0.6653 

of  n  for  9^~  1.48,  1.5,  and  1.6.  The  curve  shows  the  qua¬ 
dratic  nature  of  the  data,  and  is  useful  in  estimating  the  Hopf 
bifurcation  point,  9-\ .  The  computed  curve  is 
$7/'=  0.5 1470“ -0.00 130+  1,468,  which  yields  9''*;,=  1.468. 
This  value  is  2X10“^  lower  in  vortex  strength  than  the 
§7/-=  1  47  case  where  steady  flow  was  found.  Thus  the  Hopf 
point  identified  through  the  curve  fit  is  slightly  lower  in  vor¬ 
tex  strength  than  the  bracket  containing  the  Hopf  point  iden¬ 
tified  through  time  integration.  The  linear  term  in  the  poly¬ 
nomial  is  nearly  zero,  which  implies  that  O  departs  from 
zero  at  the  bifurcation  point  with  nearly  infinite  slope.  This  is 
in  good  agreement  with  Hopf  bifurcation  theory which 
dictates  that  the  bifurcation  parameter,  9^,  is  symmetric  with 
respect  to  a  characteristic  amplitude  of  the  bifurcated  flow, 
i.e.,  9''XCl)  =  9X-Ct).  This  symmetry  condition  holds  when 
the  linear  term  is  identically  zero. 

The  stability  of  2-D  flows  computed  with  the  noncolum- 
nar  (zero  azimuthal  vorticity)  inflow  condition  is  briefly  in¬ 
vestigated  for  Re=  1000.  The  primary  limit  point  for  this 
case  is  higher  (9"'^=  1.59048)  than  the  primary  limit  point 
for  the  columnar  inflow  condition  (9^p=  1.47435).  It  is  also 
found  that  the  region  of  nonuniqueness  for  this  case  has  in¬ 
creased  to  0.048,  compared  to  0.015  for  the  columnar  initial 
conditions.  Four  time-integration  runs  are  performed.  Two 
runs  are  performed  at  1.585^'^  and  9^=  1.585'"“,  which 
are  below  the  primary  limit  point.  The  other  two  runs  are  at 
9^=  1.595  and  9^=  1.605,  which  are  above  the  primary  limit 
point.  The  temporal  behavior  of  the  two  runs  computed  be¬ 
fore  the  primary  limit  point  is  similar  to  the  behavior  ob¬ 
served  for  the  columnar  inflow  conditions.  That  is,  both 
flows  time-asymptote  to  steady,  axisymmetric  conditions. 
Past  the  primary  limit  point,  the  9^ =  1.595  case  is  also  found 
to  approach  a  steady,  axisymmetric  flow.  It  is  noted,  how¬ 
ever,  that  much  longer  run  times  (r^  1000)  are  required  to 
approach  time-asymptotic  conditions  for  this  case.  At 
9'=  1.605,  the  initial  breakdown  bubble  transitions  to  a  3-D, 
time-periodic  flow  with  no  breakdown,  analogous  with  the 
results  for  the  columnar  inflow  condition.  Thus  the  results 
indicate  that  changing  inflow  conditions  delays  the  instability 
to  3-D  disturbances  to  a  point  slightly  past  the  primary  limit 
point.  The  difference  in  vortex  strength  between  the  primary 
limit  point  and  the  1.605  case  is  1.452X  10  Vortex 
strengths  between  these  two  cases  were  not  investigated. 

As  previously  discussed  in  Section  I,  the  location  of  the 
primary  limit  point  in  an  axisymmetric  solution  space  has 
been  associated  with  the  appearance  of  critical  flow.  The 
results  of  this  study  may  add  further  insight  into  the  role  of 
the  primary  limit  point.  The  results  show  that  the  Hopf  and 
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primary  limit  point  locations  differ  by  a  maximum  of 
5.65 X  10"^  in  Also,  for  the  case  of  noncolumnar  inflow 
conditions,  the  Hopf  point  and  the  primary  limit  point  loca¬ 
tions  differ  by  a  maximum  of  1.452X  10“^.  Thus,  the  results 
imply  a  possible  connection  between  the  loss  of  stability  to 
3-D  disturbances,  the  location  of  the  primary  limit  point,  and 
the  appearance  of  critical  flow. 

The  results  of  this  work  do  not  imply  that  the  Hopf  bi¬ 
furcation  and  the  location  of  the  primary  limit  point  are  ex¬ 
actly  coincident.  The  vortex  strengths  chosen  here  are  not 
sufficiently  close  to  the  limit  point  to  make  this  determina¬ 
tion.  Also,  subtle  differences  between  the  PAC  and  TANS 
models,  such  as  the  introduction  of  artificial  damping,  are 
unavoidable.  This  makes  it  very  difficult  to  precisely  predict 
the  onset  of  the  instability.  Furthermore,  time-integration 
proves  to  be  quite  demanding  on  computational  resources 
when  computing  solutions  near  the  Hopf  bifurcation.  This  is 
primarily  due  to  the  lengthy  times  required  for  the  instability 
to  develop,  with  this  time  increasing  dramatically  the  closer 
one  computes  to  the  actual  Hopf  point.  Instead,  a  contribu¬ 
tion  of  this  work  is  that  a  Hopf  bifurcation  point  is  bracketed 
to  lie  in  close  proximity  to  the  primary  limit  point — a  con¬ 
tribution  previously  unreported  in  the  literature.  It  is  antici¬ 
pated  that  a  more  direct  means  of  computing  the  Hopf  point 
is  in  order  for  determining  if  their  exists  a  definite  correspon¬ 
dence  between  the  flow  instability  and  the  primary  limit 
point.  Such  solvers  have  recently  been  developed  for  2-D 
flows.^"^ 


IV.  SUMMARY  AND  DISCUSSION 

Axisymmetric  solutions  computed  in  this  study  at 
Re=  1000  contain  a  region  of  nonuniqueness,  resulting  in 
primary  and  secondary  limit  points.  These  axisymmetric  so¬ 
lutions  are  treated  as  initial  conditions  to  a  3-D  time- 
integration  algorithm.  When  the  vortex  strength  is  prescribed 
to  be  less  than  that  at  the  primary  limit  point,  the  resulting 
time-asymptotic  solutions  are  found  to  be  steady  and  axi¬ 
symmetric.  In  particular,  when  three  nonunique  initial  con¬ 
ditions  are  time-integrated,  all  three  initial  conditions  time- 
asymptote  to  the  solution  corresponding  to  the  initial 
condition  of  the  upper  stable  branch.  The  transition  between 
the  initial  bubble  breakdown  solution  of  the  lower  stable 
branch  to  the  nearly  columnar  upper  stable  branch  solution  is 
investigated.  It  is  found  that  the  bubble  breakdown  flow 
transforms  temporarily  to  a  3-D  spiral  form,  which  ulti¬ 
mately  decays  to  nearly  columnar  flow.  When  the  vortex 
strength  is  prescribed  to  be  greater  than  the  primary  limit 
point,  the  resulting  time-asymptotic  solutions  are  found  to  be 
time- periodic  and  three-dimensional.  Thus  a  Hopf  bifurca¬ 
tion  is  found  to  occur  in  very  close  proximity  to  the  primary 
limit  point.  Modifying  the  inflow  boundary  condition  to  relax 
the  columnar  flow  assumption  results  in  the  Hopf  point  lying 
at  a  vortex  strength  that  is  slightly  greater  than  the  primary 
limit  point. 

Previous  results^^  by  the  authors  for  Re=250  also  show 
evidence  for  the  Hopf  bifurcation,  which  occurs  in  this  case 
near  9^ =  1.53.  However,  only  unique  axisymmetric  solutions 
exist  at  this  lower  Reynolds  number.  The  lack  of  limit  points 


preclude  any  obvious  correlation  between  the  loss  of  stability 
and  the  character  of  the  axisymmetric  solution  branch.  Two 
other  similarities  in  the  solution  behavior  are  noted  at  the 
two  Reynolds  numbers.  First,  as  9^  is  increased,  the  loss  of 
stability  occurs  before  reversed  flow  is  present.  Second,  for 
vortex  strengths  near  the  instability,  the  time-asymptotic  val¬ 
ues  of  Q  for  3-D  flows  are  higher  than  for  the  initial  2-D 
flow. 

The  role  of  the  primary  limit  point  is  further  defined  in 
this  study.  Previous  (axisymmetric)  work  had  established  a 
correspondence  between  the  primary  limit  point  and  a 
change  in  flow  criticality,  implying  that  flow  criticality  sig¬ 
nals  the  emergence  of  axisymmetric  breakdown.  This  study 
indicates  that  the  primary  limit  point  is  in  close  proximity  to 
the  vortex  strength  where  stability  to  3-D  disturbances  is 
lost,  with  further  increases  in  vortex  strength  ultimately  re¬ 
sulting  in  vortex  breakdown.  This  result  may  allow  the  dis¬ 
tinct,  time-periodic  fluctuations  associated  with  the  Hopf  bi¬ 
furcation  to  be  used  to  signal  the  onset  of  breakdown. 
However,  a  precise  statement  concerning  the  connection  be¬ 
tween  the  Hopf  bifurcation  and  the  primary  limit  point  is  not 
established  in  this  work.  This  is  due  to  the  practical  limits  of 
bracketing  a  bifurcation  point  through  time  integration,  and 
possibly  due  to  subtle  differences  in  which  artificial  damping 
is  employed  within  the  two  numerical  models. 

Finally,  it  is  of  interest  to  contrast  the  role  of  the  primary 
limit  point  for  3-D  flows  to  the  role  associated  with  axisym¬ 
metric  solutions.  For  axisymmetric  flows,  the  primary  limit 
point  is  associated  with  values  of  Q  which  are  positive. 
However,  a  slight  increase  in  above  the.  primary  limit 
point  brings  about  large  differences  in  solution  behavior,  due 
to  the  sudden  drop  of  solutions  from  the  upper  to  lower 
stable  branch.  Thus,  a  slight  increase  in  W'  from  in  2-D 
flows  results  in  a  dramatic  change  in  solution  character — 
from  flows  with  no  reversed  flow  to  flows  with  reversed 
flow.  Results  in  this  study  indicate  that  3-D  flows  do  not 
undergo  this  large  change  in  flow  structure  as  the  primary 
limit  point  is  crossed.  Instead,  flows  computed  just  past  the 
primary  limit  point  are  void  of  reversed  flow.  In  3-D  flows, 
therefore,  the  primary  limit  point  is  not  associated  with 
breakdown,  as  it  is  generally  perceived  in  axisymmetric 
flows.  For  example,  Fig.  11(a)  shows  that  the  primary  limit 
point  is  at  9^^=1.47435,  whereas  reversed  flow  is  not 
achieved  until  approximately  9^=  1.6.  In  general,  the  effect 
of  three-dimensionality  is  the  delay  of  the  formation  of  re¬ 
versed  flow,  requiring  generally  larger  values  of  vortex 
strength  to  achieve  breakdown  than  in  axisymmetric  flows. 

Further  work  is  needed  to  establish  a  precise  relationship 
between  the  Hopf  bifurcation  and  the  primary  limit  point.  In 
particular,  a  theoretical  study  of  2-D  and  3-D  vortical  flows 
in  tubes  would  help  confirm  and  may  better  define  the  nu¬ 
merically  determined  relationship  found  in  this  work.  In  ad¬ 
dition,  more  studies  are  needed  to  determine  the  effects  of 
inlet  flow  profiles  and  tube  geometry  on  the  dynamic  and 
time-asymptotic  behavior. 
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APPENDIX:  SOLUTION  PROCEDURE 


The  solution  procedure  is  a  Beam-Warming,^^  approxi- 
mate  factorization  scheme.  The  equations  are  solved  in  a 
nonconservative,  implicit  manner  using  three-point*- 
backward  time  integration.  Temporal  accuracy  is  second  or¬ 
der  for  the  inviscid  terms  and  first  order  for  the  viscous 
terms,  owing  to  the  explicit  treatment  of  the  viscous  terms.  A 
fourth-order  spatially  accurate  compact, or  Fade,  operator 
is  employed  to  approximate  derivatives.  The  compact  ap¬ 
proximation  for  the  derivative  of  a  scalar,  with  respect  to 
X,  for  example,  is  given  by 


2A.r 


^  S^ii 


where 

=  w  /  4-  [  “  w  / _  I ,  ^  +  ( tq-  + 1  ““  2  iq-  _  ^ )  /6. 

,3^  denotes  a  general  operator  while  the  subscript  refers  to  the 
coordinate  direction. 

The  general,  approximately  factored  scheme  in  delta 
form,  where  ^  -  f/'"  is  given  by 

=  (Al) 

where 
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The  and  W  matrices  at  the  unknown  time  level, 

A2  +  I ,  are  obtained  through  extrapolation.  This  is  required  for 
the  nonconservative  scheme  to  achieve  second-order  tempo¬ 
ral  accuracy.  To  insure  steady-state  consistency,'*^  the  ex¬ 
plicit  damping  coefficient,  ,  is  defined  above  such  that  the 
explicit  damping  term  is  multiplied  by  Af.  The  implicit 
damping  coefficient,  w,- ,  is  set  to  zero  for  this  work. 

The  solution  is  updated  following  f,  77,  and  f  sweeps: 

( Ar/3M''-^  •  5^]  A"  [/i 

[.y'^/+(Ar/3).^^^^^^]A«f/2=.^,(A'’f/i), 

[.^^/  +  ( A  r/3)  8  A”  f/2) . 


5^^  is  multiplied  through  Eq.  (Al)  before  solving  the  ^ 
sweep,  and  acts  on  all  of  the  right-hand  side  terms  except  the 
damping  term,  tl.  The  damping  term  is  excluded  for  consis¬ 
tency  in  comparing  results  to  the  central-difference  PAC 
scheme.  The  traditional  method  of  sweeping  in  the  77  and  i 
directions  is  modified^^  to  allow  for  the  presence  of  a  multi¬ 
block  grid  structure.  Spatial  derivatives  in  f3,  grid  metrics, 
and  viscous  terms,  D,  are  computed  implicitly  at  time-level 
n  to  fourth-order  accuracy.  This  is  accomplished  in  an  effi¬ 
cient  manner  using  a  lower/upper  decomposition  technique, 
since  the  resulting  set  of  equations  represent  a  constant,  tridi¬ 
agonal  matrix-inversion  problem.  At  boundaries  si,  s2  and 
s3,  second-order  accurate,  three-point  forward/backward  dif¬ 
ferences  are  used.  Derivative  terms  normal  to  the  grid  cuts 
are  approximated  with  fourth-order  central  differences. 

The  numerical  procedure  was  partially  checked  by  solv¬ 
ing  for  a  uniform  freestream  with  Dirichlet  (no  correction) 
boundary  conditions  at  si,  s2  and  s3.  The  resulting  temporal 
corrections,  A"f/,  were  found  to  be  machine  zero  every¬ 
where  for  100  iterations  over  a  wide  range  of  time  steps.  It  is 
of  interest  to  note  that  employing  a  fourth-order  accurate 
boundary  stencil^^  along  the  physical  boundaries  si,  s2  and 
s3  of  the  form 

1  +  3m^2“  (  17Ki  +  9(«2  +  “3)“^^4) 

resulted  in  freestream  conservation  errors. 

Validation  of  the  TANS  model“®  proceeded  along  two 
lines.  First,  a  variety  of  model  problems  were  solved,  includ¬ 
ing  the  incompressible  and  compressible  flows  over  a  flat 
plate,  unsteady  Couette  flow  and  unsteady  heat  conduction. 
The  model  output  and  the  known  solutions  matched  to  a  high 
degree  of  accuracy  in  each  case.  Second,  the  TANS  and  PAC 
models  were  cross-validated  by  comparing  computed  solu¬ 
tions  obtained  under  flow  conditions  which  lead  to  steady 
and  axisymmetric  flow.  The  resulting  solutions  were  found 
to  agree  to  a  high  degree  of  accuracy. 
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Force  and  moment  experiments  were  conducted  for  a  65-deg  delta  wing  undergoing  ramp-and-hold 
and  harmonic  rolling  motions.  This  extensive  set  of  experiments  isolated  critical-state  responses.  Motions 
between  critical  states  and  motions  crossing  critical  states  were  included.  Roll-angle  amplitude  and  the 
roll  rate  were  varied.  The  total  angle  of  attack  and  Mach  number  were  held  constant  at  30  deg  and  0.3, 
respectively.  The  amount  of  time  required  for  the  rolling  moment  coefficient  to  reach  its  steady-state  value 
after  the  end  of  a  ramp  was  quantified  for  numerous  ramp  motions.  This  relaxation  time  was  a  signifi¬ 
cantly  large  value  for  many  of  the  motions,  especially  when  the  5-deg  critical  state  was  crossed.  Motion 
history  effects  were  determined  for  motions  that  crossed  critical  states.  The  effect  of  roll  rate  on  critical- 
state  transients  proved  to  be  insignificant. 


Nomenclature 

b  =  trailing-edge  wing  span 
Cl  -  rolling  moment 

Cl  ~  Cl  Cl 

'dyn  ‘  '(Otic 

c  -  root  chord 
c  =  mean  aerodynamic  chord 
k  =  reduced  frequency,  ojbHU^ 
q  =  freestream  dynamic  pressure 
S  =  model  planform  area 
t  -  time 

=  nondimensional  time,  lU^t/b 
f/co  -  freestream  velocity 

cr  =  total  angle  of  attack,  body-axis  inclination  with 
respect  to 

-  body-axis  roll  angle  ^ 

0  =  first  derivative  of  </>  with  respect  to  time 

=  second  derivative  of  <f)  with  respect  to  time 
=  nondimensional  roll  rate, 

4>i  =  shaft  rotation  angle  at  the  base  of  the  sting 

02  =  shaft  rotation  angle  at  the  center  of  the  force  balance 

CO  =  circular  frequency,  rad/s 

Introduction 

N  expansion  of  the  flight  envelope  has  been  and  remains 
a  very  important  goal  for  designers  of  high-performance 
air  vehicles.  Such  high-performance  aircraft  are  required  to  fly 
faster  and  at  higher  angles  of  attack.  Because  of  these  demands 
they  often  employ  highly  swept  wing  surfaces.  When  flying  at 
higher  angles  of  attack,  controllability  has  become  a  large 
problem  for  these  aircraft.  Intense  changes  in  the  flowfield  are 
encountered  that  result  in  dramatic  nonlinearities  in  an  air¬ 
craft’s  aerodynamic  behavior.  Under  these  conditions  the  aero¬ 
dynamics  are  time-  and  motion-history-dependent;  therefore, 
conventional  quasisteady  modeling  techniques  are  not  ade¬ 
quate. 

This  paper  will  discuss  some  implications  of  experimental 
force  and  moment  results  for  a  65-deg  delta  wing  driven 
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through  roll  maneuvers.  This  delta  wing  has  been  the  focus  of 
an  in-depth  study  over  the  past  several  years.*" ’’  Previous 
analyses  of  the  delta  wing  have  shown  severe  nonlinear  aero¬ 
dynamic  behavior  and  the  aerodynamics  have  been  shown  to 
be  dependent  on  the  precise  motion  history.  As  with  other 
highly  swept  planforms,  vortex  breakdown  was  found  to  be  a 
large  contributor  to  the  aerodynamic  behavior  of  the  65-deg 
delta  wing.  The  axial  progression  of  breakdown  because  of  a 
time-dependent  variation  in  wing  orientation  was  found  to  be 
quite  slow. 

The  dynamic  force  and  moment  responses  possessed  large 
time  lags.  Two  possible  causes  for  these  time  lags  are  1)  finite 
response  times  caused  primarily  by  the  response  of  vortex 
breakdown  and  2)  transients  resulting  from  one  or  more  crit¬ 
ical-state  encounters.  A  critical  state  exists  when  a  disconti¬ 
nuity  in  a  static  force  or  moment  curve  (or  its  derivative)  oc¬ 
curs  and  its  crossing  causes  a  discrete  change  in  the 
equilibriuiti  response.  Thus,  time  lags  can  result  without  en¬ 
countering  a  critical  state,  but  if  a  critical  state  is  encountered, 
there  is  an  additional  amount  of  time  that  is  required  for  the 
transition  from  one  flow  state  to  another.  The  time  scales  for 
these  critical-state  transients  had  not  been  established.  For 
modeling  purposes  and  possible  simplifications,  it  would  be 
very  useful  to  know  if  the  effects  from  one  or  more  critical- 
state  encounters  can  be  neglected. 

The  force  and  moment  experiments  of  this  study  were  de¬ 
signed  to  isolate  the  effects  induced  by  critical-state  encoun¬ 
ters.  Measurements  were  obtained  for  ramp-and-hold  motions 
and  harmonic  motions  in  roll.  Comprehensive  static  and  dy¬ 
namic  force  and  moment  measurements  were  obtained.  An  ex¬ 
amination  of  these  results  will  be  the  focus  of  this  paper,  with 
the  goal  being  to  increase  the  knowledge  of  how  to  better 
model  the  transients  induced  from  critical-state  encounters. 

Discussion 

Experimental  Equipment 

The  model  used  for  the  experiments  was  a  65-deg  delta  wing 
with  sharp  leading  and  trailing  edges.  A  schematic  of  the 
model  is  shown  in  Fig.  1.  Both  the  leading  and  trailing  edges 
were  symmetrically  beveled  with  an  included  angle  normal  to 
the  leading  edge  of  20  deg.  The  wing  was  constructed  with  a 
multilayer  graphite  composite  skin  and  a  foam  core,  and  the 
fuselage  was  constructed  purely  of  graphite  composite.  The 
wing  thickness-to-root-chord  ratio  is  1.53%. 

The  delta  wing  was  mounted  in  the  Subsonic  Aerodynamic 
Research  Laboratory  (SARL)  7  ft  X  10  ft  (2.1  m  X  3.0  m) 
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Sharp  Leading  and 
Trailing  Edges 


Symmetric  Bevel 
(20°  Included  Angle) 


Area=:1.94 

(0.180  m2) 


.i/c=\,53% 


Fig.  1  Delta  wing  model. 


Fig.  2  Roll  rig  apparatus  in  SARL  wind  tunnel. 


wind  tunnel  via  a  roll  rig.  The  wing  was  mounted  on  a  five- 
component  strain  gauge  balance  that  was  supported,  by  a  drive 
shaft,  through  a  hollow  sting.  The  dynamic  roll  rig,  installed 
in  the  SARL  tunnel,  is  shown  in  Fig.  2.  The  apparatus,  devel¬ 
oped  by  the  Canadian  National  Aeronautical  Establishment 
(NAE),  was  designed  for  large-amplitude  and  high-rate  mo¬ 
tions.  A  more  detailed  description  of  both  the  model  and  roll 
rig  design  can  be  found  in  Ref.  12. 

Data  Acquisition  and  Reduction 

The  data  presented  here  were  obtained  from  a  series  of  ex¬ 
periments  designed  and  conducted  both  by  personnel  from 
Wright  Laboratory  and  members  of  the  Canadian  Institute  for 
Aerospace  Research  (lAR)  in  1994,  The  delta  wing  was  driven 
through  both  ramp-and-hold  maneuvers  and  harmonic  oscil¬ 
lations  in  roll.  During  these  experiments,  digitally  sampled 
time  histories  of  both  the  model  motion  and  five-component 
force  and  moment  responses  were  obtained. 

For  all  roll  motions  presented  here  the  total  angle  of  attack 
cr  was  fixed  at  30  deg  and  the  Mach  number  was  0.3.  Data 
were  acquired  with  a  constant  sampling  interval  for  a  given 
motion  and  for  several  (usually  20)  cycles  of  the  motion.  Five- 
component  force  and  moment  data  were  ensemble-averaged  to 
yield  1024  data  points  for  one  cycle.  The  five  components 
produced  were  normal  force,  pitching  moment,  side  force, 
yawing  moment,  and  rolling  moment.  The  shaft  rotation  angle 
at  the  base  of  the  sting  was  also  measured  throughout  the  same 
sampling  time.  For  each  motion  data  were  acquired  both  with 
the  tunnel  operating  at  the  desired  testing  speed  and  with  the 
wind  off. 

Ramp  Motions 

The  constant  ramp  roll  rate  was  varied  from  0.25  rad/s  (14 
deg/s)  through  7  rad/s  (400  deg/s)  for  the  various  ramp  mo¬ 
tions.  A  ramp-up  cycle  consists  of  a  hold  at  the  ramp- 
up,  a  0.5-s  hold  at  and  a  return  to  The  data,  how¬ 
ever,  were  acquired  only  for  a  small  amount  of  time  before 


the  ramp  up  starts,  through  the  ramp  up,  and  for  the  0.5-s  hold 
after  the  end  of  the  ramp.  The  part  of  the  cycle  spent 
returning  to  4>m\n  ^nd  an  amount  of  time  at  before  starting 
the  next  cycle  are  not  included  in  the  data  acquisition  sampling 
period.  Data  were  acquired  over  a  sufficient  number  of  cycles 
so  that  the  ensemble-averaged  data  did  not  change  with  the 
acquisition  of  data  over  an  additional  cycle. 

The  first  step  taken  in  the  data  reduction  process  was  to 
convert^  the  raw  output  voltages  from  the  force  balance  to  non- 
dimensional  force  and  moment  coefficients,  based  on  the 
model  geometry.  The  remaining  part  of  the  data  reduction  was 
much  more  involved.  For  many  applications,  a  simple  sub¬ 
traction  of  the  wind-off  tare  data  from  the  wind-on  data  would 
be  an  adequate  method  of  removing  the  inertial  effects.  How¬ 
ever,  for  the  ramp  motions  at  which  high  angular  accelerations 
were  commanded  at  the  beginning  and  end  of  the  ramp,  it  was 
determined  that  such  a  subtraction  was  not  appropriate. 

For  a  given  ramp-and-hold  motion,  the  difference  in  the  roll 
angle  between  the  wind-on  and  wind-off  data  was  very  small. 
Slight  differences  between  the  measured  and  commanded  roll 
angles,  however,  resulted  in  severely  different  accelerations 
commanded  by  the  servo-system.  Because  of  this,  significant 
oscillations  were  present  in  both  the  wind-on  and  wind-off 
rolling  moment  coefficients,  particularly  following  the  com¬ 
manded  accelerations.  The  oscillations  between  the  wind-on 
and  wind-off  runs  were  not  always  in-phase  because  the  small 
roll-angle  errors  are  influenced  by  the  airloads  acting  on  the 
model.  Therefore,  a  subtraction  of  the  two  rolling-moment 
time  histories  often  increased  the  amplitude  of  the  oscillations. 

The  inertial  effects  were  determined  by  an  alternate  method 
and  they  were  subtracted  from  the  wind-on  force  and  moment 
coefficients,  yielding  a  more  accurate  result,  with  significantly 
smaller  oscillations.  Knowing  the  roll  angle,  the  rolling  mo¬ 
ment  contribution  caused  just  by  the  inertia  can  be  calculated 
and  subtracted  from  the  wind-on  rolling  moment  resulting 
from  both  the  inertial  and  aerodynamic  effects. 

To  determine  the  inertial  effects  the  roll  motion  history  <^,(r) 
from  the  tare  runs  was  used  to  calculate  <5^,(0.  ^^(t)  was  cal¬ 
culated  by  taking  a  second-order  central  difference  in  time. 
The  initial  and  final  values  of  ^,(r)  were  calculated  using  a 
second-order  one-sided  difference.  Because  of  the  noise  pro¬ 
duced  from  differentiating  twice,  frequencies  in  the  roll  accel¬ 
eration  time  histories  above  100  Hz  were  numerically  filtered 
out.  Knowing  0,(/)  and  C/(/)  for  a  tare  run,  an  estimate  of  the 
inertia  I  could  be  determined  assuming  a  simple  one-degree- 
of-freedom  system  in  roll  where 


Ci(t)  =  (f/qSb)^i(t)  +  AT  (1) 

K  is  the  value  of  Ci  before  and  after  the  ramp  when  0,(r)  =  0. 
This  value  was  near  zero,  but  is  included  to  account  for  any 
small  biases  in  the  signal.  Using  a  least-squares  linear  curve 
fit,  values  of  both  I  and  K  were  calculated,  since  q  from  the 
wind-on  tests  S  and  b  are  all  known.  Also,  the  correlation 
coefficient  could  be  calculated  to  assess  the  validity  of  the 
linear  correlation  assumption.  Correlation  coefficients  from 
0.60  to  0.85  were  obtained.  Plots  of  C,  vs  0,  showed  noticeable 
loops  around  the  linear  curve  fit.  Since  the  rotation  angle  was 
sensed  at  the  end  of  the  actuator  shaft,  it  was  desired  to  try  to 
determine  the  shaft  rotation  angle  at  the  same  location  where 
the  rolling  moment  was  measured,  at  the  balance  center.  It  was 
thought  that  a  difference  in  the  rotation  angle  between  the  base 
of  the  shaft  and  the  balance  center,  because  of  torsional  de¬ 
flections  of  the  drive  shaft,  might  be  the  cause  of  the  observed 
loops. 

To  determine  the  shaft  rotation  angle  at  the  balance  center 
4>2  the  entire  shaft  from  its  base  to  the  center  of  the  balance 
was  modeled  as  a  cylinder  of  uniform  torsional  stiffness.  This 
system  was  modeled  as  a  second-order  system  in  terms  of  the 
difference  in  shaft  rotation  angle  A<^,  where  A0 
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Therefore,  the  system  was  described  by  the  following  equa¬ 
tion: 

A<ji(0  +  2^w„A<^(f)  +  <olA(f>U)  =  0  (2) 

(  and  oj„  are  the  damping  ratio  and  natural  frequency  of  the 
drive  shaft  torsional  degree  of  freedom,  respectively.  A<f)(t) 
was  calculated  using  the  indicial  form  of  the  Duhamel  integral: 

Am  =  Ut)  -  =  A(t)<t>m  +  1^0-  »<A(t)  dr  (3) 

where  =  —g'^Tcos  cot  +  (/T/£j)sin  cot]  with  n  =  and 
CO  =  a>„Vl  “  Since  <^,(r)  was  known,  <^2(0  was  then  de¬ 
termined.  (pi(t)  was  found  to  be  essentially  identical  to  (f>2(t). 

Fixed  values  of  {  and  co„  were  found  by  trying  to  minimize 
the  difference  between  the  measured  Cj(t)  and  the  Q(t)  cal¬ 
culated  using  ^2(0  in  Eq.  (1).  Previous  work"  had  shown, 
however,  that  the  damping  force  was  not  proportional  to 
but  rather  it  was  essentially  a  constant  force  opposing  the  di¬ 
rection  of  motion.  Therefore,  f  was  set  to  be  inversely  related 
to  I  I  to  provide  a  constant  damping  force.  This  provided  a 
good  match  between  the  measured  and  calculated  values  of 
Ci(t)  as  desired. 

Using  the  values  of  I/qSb  and  K  determined  by  including 
the  tare  values  of  in  Eq.  (1),  the  necessary  information 
was  found  to  subtract  the  inertial  effects  from  the  wind-on 
rolling  moment  data.  For  ramp  motions,  the  values  of  C/(r) 
discussed  in  the  remainder  of  this  paper  were  calculated  by 
taking  the  wind-on  measured  values  of  C/(t)  and  subtracting 
C{(t)  calculated  from  Eq,  (1),  using  the  wind-on  measured  val¬ 
ues  for  (f>i(t). 

Harmonic  Motions 

The  frequency  was  varied  from  1 . 1  through  1 1  Hz  for  the 
harmonic  motions.  For  these  motions  one  sampling  cycle  is 
one  complete  harmonic  oscillation.  In  these  cases,  force  and 
moment  data  were  obtained  by  subtracting  the  wind-off  data 
from  the  wind-on  data.  Obvious  oscillations  were  present  in 
the  resulting  force  and  moment  time  histories.  To  better 
smooth  the  data,  a  fast  Fourier  transform  was  performed  and 
a  limited  number  of  harmonics  were  kept.  The  number  of  har¬ 
monics  retained  varied  with  forcing  frequency  as  shown  in 
Table  1.  Determining  the  number  of  harmonics  to  retain  was 
established  by  choosing  a  number  small  enough  to  provide 
sufficient  smoothing  while  retaining  the  overall  shape  of  the 
time  histories. 

Observations 

Previous  findings^  revealed  the  existence  of  a  critical  state 
for  this  delta  wing  at  approximately  =  5  deg.  This  can  be 
seen  in  the  static  rolling  moment  coefficient  vs  0,  shown  in 
Fig.  3.  In  the  region  of  =  5  deg  a  definite  discontinuity  in 
the  rolling  moment  coefficient  is  observed.  A  discontinuity  in 


Table  1  Number  of  harmonics  retained  for  harmonic 
motions  of  varying  forcing  frequency 


Forcing 

frequency, 

Hz 

Reduced 

frequency 

Number  of 
harmonics 
kept 

Resulting 

cutoff 

frequency, 

Hz 

1.1 

0.010 

8 

8.8 

2.2 

0.021 

8 

17.6 

3.3 

0.031 

7 

23.1 

4.4 

0.042 

6 

26.4 

5.5 

0.052 

5 

27.5 

7.7 

0.073 

4 

30.8 

8.8 

0.084 

4 

35.2 

11 

0.105 

3 

33.0 

the  pitching  moment  coefficient  was  also  observed  around  this 
angle.  These  discontinuities  indicate  the  existence  of  a  sub- 
critical  bifurcation.  Critical  states  around  8  and  1 1 .3  deg  were 
also  proposed  in  previous  work.  The  8 -deg  critical  state  is 
likely  to  be  a  Hopf  bifurcation  because  of  the  large  rms  fluc¬ 
tuations.^  The  increase  in  scatter  for  ~8  deg  <  </>  <  ~12  deg, 
shown  in  Fig.  3,  indicates  that  time-averaging  failed  to  average 
out  low-frequency  fluctuations.  The  11.3-deg  critical  state  ap¬ 
pears  to  be  a  supercritical  bifurcation.  Also,  the  unsteadiness 
in  Cl  is  greatly  reduced  for  <^  >  —12  deg. 

Ramp  Motions 

The  experiments  were  designed  to  isolate  the  effects  of  crit¬ 
ical-state  encounters  from  viscous  time  lags  caused  primarily 
by  the  slow  response  of  vortex  breakdown  position.  Four  dif¬ 
ferent  ramp  motions  are  shown  in  Fig.  4.  These  ramp  motions 
either  stay  below  the  critical  state  around  <^  =  5  deg,  cross  the 
critical  state,  or  stay  above  the  critical  state. 

Figure  5  illustrates  CXO  for  a  ramp  motion  in  roll  from  7  to 
4  deg  at  —1  rad/s  as  compared  to  the  quasisteady  C/(0  over 
this  range  in  roll  angle.  This  figure  reveals  an  obvious  differ¬ 
ence  between  the  measured  and  quasisteady  response.  Also,  it 
is  interesting  that  at  ramp  onset  the  dynamic  response  is  in¬ 
creasing,  while  the  reverse  is  true  of  the  quasisteady  behavior. 
Note  that  C/(0  does  not  approach  the  quasisteady  value  until 
long  after  the  ramp  has  ended.  The  amount  of  time  required 
for  the  rolling  moment  coefficient  after  the  end  of  a  ramp  (<^ 
0)  to  reach  its  quasisteady  value  was  quantified  for  many 
motions. 

It  was  desired  to  isolate  the  transient  effects  resulting  from 
the  presence  of  the  leading-edge  vortices.  Previous  findings’^ 


Fig.  4  Sample  (1  rad/s)  ramp  motions  about  the  5-deg  critical 
state. 
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Fig.  5  Static  and  dynamic  Q  for  a  (-1  rad/s)  ramp  through  the 
5-deg  critical  state. 


Fig.  6  Determination  of  relaxation  times  for  ramp  motions 
through  the  5-deg  critical  state  (<^nd  *  0.003). 


revealed  that  the  vortical  component  of  the  rolling  moment 
had  a  net  destabilizing  effect  for  roll-angle  variations.  The 
presence  of  vortex  breakdown  and  its  relatively  slow  response 
is  considered  to  be  the  dominant  cause  of  time  lags  observed 
for  dynamic  maneuvers  (not  counting  critical-state  effects). 
The  potential  flow  component  of  the  rolling  moment  reacts  at 
the  freestream  convection  speed,  which  for  these  rates  can  be 
considered  virtually  instantaneous.  Therefore,  the  potential 
flow  part  of  the  rolling  moment  coefficient  was  subtracted 
from  the  total  rolling  moment  coefficient.  The  potential  con¬ 
tribution  was  calculated  using  QUADPAN,  a  panel  code  de¬ 
veloped  by  Lockheed-California  Company.*"^ 

Figure  6  shows  the  vortical  contribution  to  the  rolling  mo¬ 
ment  coefficient,  the  total  rolling  moment  minus  the  potential 
contribution.  This  figure  illustrates  how  relaxation  times  were 
calculated.  The  first,  second,  and  third  rolling  moment  coeffi¬ 
cients  correspond  to  ramp  motions  from  “4  deg  to  4,  6,  and 
7  deg,  respectively,  each  for  a  roll  rate  of  1  rad/s,  as  shown  in 
Fig.  4.  For  each  of  these  three  ramp  motions  the  vortical  roll¬ 
ing  moment  responses  were  time-shifted  so  that  a  nondimen- 
sional  time  of  0.0  corresponded  to  the  end  of  each  of  the  ramp 
motions.  Also,  the  second  vortical  rolling  moment  coefficient 
(AC/)2,  shown  in  Fig.  6,  has  the  response  of  the  first  motion 
subtracted.  Likewise,  the  third  response  had  the  response  from 
the  second  motion  subtracted. 

The  relaxation  times  that  were  quantified,  as  shown  in  Fig. 
6,  were  actually  the  time  after  the  end  of  the  ramp  until  the 
rolling  moment  coefficient  reached  63%  of  its  value  at  the  end 


of  the  ramp-and-hold.  This  value  of  63%  was  chosen  because, 
for  an  exponential  decaying  function  C/(r),  where  AC,(0  = 
Cwtl  “  with  c,(f  ->  00)  =  Ci^,  and  C,(t  =  a)  = 

0.632C/^^.  a  must  have  the  same  units  as  t  and  it  is  called  the 
time  constant.  The  relaxation  times,  according  to  this  defini¬ 
tion,  are  indicated  by  vertical  lines.  It  should  be  noted,  as  ev¬ 
idenced  by  Fig.  6,  that  these  calculations  of  the  relaxation 
times  are  conservative  because  they  did  not  take  the  time- 
averaged  static  values  (denoted  by  solid  symbols)  to  be  the 
necessary  equilibrium  value.  The  second  ramp  motion,  which 
crosses  the  critical  state,  illustrates  that  even  at  the  end  of  the 
sampling  period  the  vortical  rolling  moment  contribution  has 
not  yet  reached  the  static  value. 

These  (63%)  relaxation  times  were  quantified  for  several 
ramp  motions  with  positive  and  negative  roll  rates  about  var¬ 
ious  critical  states.  The  nondimensional  relaxation  times 
as  well  as  the  AC/^  (from  the  end  of  the  ramp  to  the  assumed 
steady-state  value)  are  included  in  Table  2. 

Table  2  reveals  the  relaxation  time  constant  about  the  5-deg 
critical  state  (r*^  70)  to  be  quite  large.  The  8-deg  critical 

state  yielded  a  lower  value  (f ^  50).  Calculations  of  the 
relaxation  time  constants  for  the  11.3-deg  critical  state,  al¬ 
though  not  included  in  the  table,  found  10.  The  AC/^ 

values  used  to  calculate  the  relaxation  time  constants  for  the 
11.3-deg  critical  state,  however,  were  very  small.  Therefore, 
these  relaxation  time  constants  were  less  meaningful.  This 
might  be  expected  from  Fig.  3,  which  shows  that  for  roll-angle 
variations  around  11.3  deg,  the  difference  in  the  static  rolling- 
moment  coefficient  is  small. 

The  relaxation  time  constants  calculated  for  ramp  motions 
that  did  not  involve  crossing  a  critical  state  yielded  values  of 
^63%  ^  30.  For  these  motions,  where  vortex  breakdown  existed 
over  the  wing,  the  time  constant  reflects  the  relaxation  time 
resulting  from  the  response  of  vortex  breakdown. 

The  range  of  ramp  motions  examined  allowed  history  effects 
across  critical  states  to  be  assessed.  As  shown  in  Fig.  4,  ramp 
motions  were  conducted  for  cf>  between  —4  and  6  deg,  between 
-4  and  7  deg,  and  between  6  and  7  deg.  The  rolling-moment 
data  obtained  for  the  ramp  from  6  to  7  deg  were  compared  to 
the  rolling  moment  data  for  the  ramp  from  “4  to  7  deg  with 
the  rolling- moment  data  for  the  ramp  from  —4  to  6  deg  sub¬ 
tracted.  The  result  of  this  subtraction  is  denoted  A  =  6  to  7  deg. 
This  comparison  is  shown  in  Fig.  7a.  Note  again  that  the  po¬ 
tential  flow  contribution  to  the  rolling  moment  was  subtracted, 
leaving  just  the  vortical  contribution.  This  plot  reveals  no  sig¬ 
nificant  history  effects  on  the  rolling  moment  coefficient  in  the 
</>  range  from  6  to  7  deg  after  crossing  the  5-deg  critical  state. 
Table  2  also  illustrates  that  the  time  constants  30)  de¬ 

termined  for  A  =  6  to  7  deg  were  approximately  equal  to  the 
time  constant  determined  for  a  direct  ramp  motion  from  6  to  7 
deg. 

The  same  type  of  procedure  was  used  to  examine  the  history 
effect  on  the  subsequent  response  after  passing  through  the  8- 
deg  critical  state.  Figure  7b  illustrates  the  comparison  between 
the  vortical  rolling  moment  coefficient  for  a  ramp  from  9  to 
10  deg  with  a  ramp  from  6  to  10  deg,  subtracting  the  hi.story 
of  the  ramp  from  6  to  9  deg.  This  comparison  illustrates  a 
noticeable  difference  between  the  two  results.  Therefore,  in¬ 
cluding  the  flow  history  from  6  to  9  deg,  which  crosses  the  8- 
deg  critical  state,  produced  a  distinctly  different  transient  re¬ 
sponse  over  the  9-  to  10-deg  range.  The  causes  for  and 
conditions  under  which  critical-state  history  effects  become 
important  are  not  yet  understood.  This  will  be  a  subject  of 
future  investigations. 

The  effect  of  roll  rate  on  the  critical-state  transients  was 
investigated  using  the  data  obtained  with  different  roll  rates 
over  the  same  roll-angle  range.  Figure  8,  which  contains  the 
vortical  contribution  to  the  rolling  moment  coefficient  for  three 
different  roll  rates,  reveals  that  the  roll  rate  has  no  significant 
effect  on  the  transients  for  the  5-deg  critical  state.  Again,  Table 
2  illustrates  that  the  time  constants  for  each  roll  rate  were  all 
approximately  ^  70. 
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Table  2  Relaxation  times  to  63  %  of  steady>state  values  for  various  ramp  motions 


Critical 

state, 

deg 

Rate 

rad/s, 

ND 

Up/down 

Precritical  state 

Crossing  critical  state 

Postcritical  state 

5 

—4  to  4  deg 

—4  to  6  deg 

—4  to  7  deg 

A  =  4  to  6  deg® 

A  =  6  to  7  deg 

1  (0.003) 

Up 

36 

0,010 

62  0.014 

30  0.006 

4  (0.012) 

Up 

33 

0,014 

68  0.014 

40  0.007 

7  (0.020) 

Up 

33 

0,014 

73  0.006 

47  0.007 

7  to  6  deg 

7  to  4  deg 

7  to  —4  deg 

A  =  6  to  4  deg 

A  =  4  to  —4  deg 

1  (0.003) 

Down 

24 

-0.006 

61  -0.014 

78  -0.010 

8 

6  to  7  deg 

6  to  9  deg 

6  to  10  deg 

A  =  7  to  9  deg 

A  =  9  to  10  deg 

1  (0.003) 

Up 

33 

0.006 

42  0.009 

5  0.004 

10  to  9  deg 

10  to  7  deg 

10  to  6  deg 

A  =  9  to  7  deg 

A  -  7  to  6  deg 

1  (0.003) 

Down 

14 

-0.002 

53  -0.008 

34  -0.005 

*A  5=  4  to  6  deg  denotes  moment  response  for  4>  from  —4  to  4  deg  was  subtracted  from  moment  response  for  <j)  from  —4  to  6  deg. 


End  of  Ramp  Motion 


Fig.  7  Response  of  the  vortical  contribution  to  Q  with  and  with¬ 
out  the  prior  motion  history  through  the  a)  5-  and  b)  8-deg  critical 
state  included  “  0.003). 


Harmonic  Motions 

The  rolling-moment  responses  for  numerous  harmonic  roll 
oscillations  were  obtained.  Figure  9  contains  rolling-moment- 
coefficient  responses  for  5-deg  amplitude  harmonic  oscillations 
in  roll.  These  oscillations,  which  vary  in  frequency  from  1.1 
to  11  Hz  (0.010  <  /:  <  0.105),  are  centered  about  <}>  =  3  deg. 
The  solid  symbols  in  Fig.  9  show  the  static  results  and  the 
dashed  curve  shows  the  static  potential  flow  contribution  to 
Cl.  These  C,  loops  run  in  a  counterclockwise  sense,  therefore, 


Fig.  8  Roll-rate  effect  on  the  transient  C/  response  through  the 
5-deg  critical  state. 


at  a  given  C;  on  the  upstroke  is  less  than  C/  on  the  down- 
stroke. 

For  the  oscillation  of  lowest  frequency,  k  =  0.010,  C/  at¬ 
tempts  to  follow  the  curvature  of  the  static  data.  As  the  fre¬ 
quency  is  increased,  however,  the  deviation  from  the  static  data 
becomes  more  apparent.  For  all  but  the  lowest  frequency  the 
rolling  moment  response  no  longer  even  remotely  resembles 
the  curvature  of  the  static  data.  For  each  oscillation  shown  in 
Fig.  9,  Cl  for  and  <^>max  were  obtained,  A  line  (dash-double¬ 
dot)  is  shown  connecting  these  two  values  for  each  oscillation. 
Even  at  the  lowest  frequency,  C/  values  at  the  roll  angle  limits 
(0  0)  do  not  reach  their  static  values.  With  increasing  fre¬ 

quency  the  slope  of  the  line  connecting  the  roll  angle  limits 
moves  closer  to  the  slope  of  the  potential  flow  contribution.  It 
would  be  expected  that  for  a  harmonic  oscillation  of  very  large 
frequency,  the  flow  would  respond  only  to  the  virtually  in¬ 
stantaneously  reacting  flow  component.  Furthermore,  the  am¬ 
plitude  of  the  rolling-moment  response  loops  decreases  as  fre¬ 
quency  increases.  Figure  9  also  illustrates  that  on  the  upstroke, 
in  the  range  of  the  offset  angle,  (f)  =  3  deg,  the  value  of  C,  is 
virtually  independent  of  roll  rate.  However,  on  the  downstroke, 
after  crossing  the  5-deg  critical  state,  the  roll  rate  has  a  sig¬ 
nificant  effect  on  the  rolling-moment  coefficient. 

Using  the  nonlinear  indicial  response  model  (NIR)  and  the 
superposition  integral  (presented  in  detail  in  previous  pa- 
pers*^’‘‘^),  the  classic  Taylor-series  expansion  for  the  aerody¬ 
namic  rjesponse  to  a  motion  input  can  be  derived.  When  these 
results  are  specialized  to  the  case  of  forced  harmonic  motion, 
relationships  are  established  between  specific  stability  deriva¬ 
tives  (both  linear  and  nonlinear)  and  the  variation  of  the  aero- 
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Decreasing  ^ 


Staady  Vaiuat 


Fig.  9  Cl  response  to  S-deg-ampIitude  harmonic  motions  cen¬ 
tered  at  “  3  deg. 


dynamic  responses  at  harmonics  of  the  forcing  frequency.  That 
iSt  for  the  steady-state  response  to  a  harmonic  motion  input;  e.g., 

<^(0  =  </)o  +  cos(kt)  (4) 

the  rolling  moment  response  can  be  shown  to  be  of  the  form 

Ci^  =  C/g  +  F|  cos(kt)  +  F2COs(2kt)  +  cosQkt)  H - 

+  G,  sin(kt)  +  G2  sin(2A:r)  +  G3  sinOkt)  +  •  •  •  (5) 

where  each  of  the  in-phase  coefficients  F„  is  a  summation  of 
the  contributions  f„j  from  particular  stability  parameters  that 
are  designated  by  the  subscript  j\  Therefore, 

>1 


Similarly,  out-of-phase  terms  G„  also  consist  of  contributions 
g„j.  The  individual  contributions  can  be  factored  into  the  form 

A;  A(^)[a2A:"  +  aX  +  •••] 

(6) 

8nJ  =  +  tX  +  +  *  *  ’] 

Note  that  the  frequency  effect  appears  as  a  power  series  in  L 
Also,  the  coefficients  of  these  terms,  (a„  Z?,),  correspond  to  the 
relevant  stability  derivative.  Furthermore,  the  higher-order 
terms  (in  k)  arise  from  repeated  differentiation  of  the  motion 
parameters.’^  Differentiation  of  sin(nkt)  gives  nk  cos(/2k0-  That 
is,  the  series  proceeds  in  the  order  bi(nk)  4-  a2(nkf  +  •  •  • . 
Thus,  if 


k  >  \l{na2fbi) 


the  series  diverges. 

It  must  also  be  noted  that  this  analysis  is  valid  under  the 
following  conditions: 

1)  The  motion  is  analytic  (which  is  true  for  forced  oscillation 
tests). 

2)  The  motion  does  not  cross  a  critical  state. 

3)  Time  constants  for  the  indicial  response  transients  are  less 
than  link,  where  nk  is  the  frequency  of  the  affected  harmonic. 

The  convergence  properties  of  the  series  representation 
(when  applied  to  harmonic  motion)  can  be  used  to  provide  an 
independent  check  of  indicial  response  time  scales.  From 
ramp-and-hold  motion  data  (see  Table  2),  indicial  response 
time  constants  of  about  30  are  expected  in  the  roll  angle  range 


Table  3  Regression  results 


J 

Parameter 

Harmonic, 

n 

najhx 

^crii 

1 

Linear,  0 

1 

28.4 

0.035 

2 

1 

19.5 

0.051 

3 

<l>4><t> 

1 

16.0 

0.063 

4 

2 

24.8 

0.040 

a)  Roll  Angle  [deg] 


b)  Roll  Angle  [deg] 


Fig.  10  C/j  for  harmonic  oscillations  centered  at  “  0  deg  for 
A:  =  a)  0.02  and  b)  0.04. 

between  the  critical  states  at  ±5  deg.  Thus,  for  harmonic  mo¬ 
tions  in  this  range,  the  series  approximation  is  expected  to 
converge  for  reduced  frequencies  below  about  0.033.  There¬ 
fore,  harmonic  motions  (centered  about  0-deg  roll  angle  and 
with  amplitudes  of  2,  3,  and  4  deg)  at  the  two  lowest  frequen¬ 
cies,  k  =  0.02  and  0.04,  were  analyzed.  A  stepwise  regression 
analysis’  was  used  to  identify  the  most  significant  terms  in  the 
series  expansion.  The  correlation  coefficient  for  the  regression 
was  0.886,  and  comparisons  between  experimental  and  pre¬ 
dicted  are  shown  in  Figs.  10a  and  10b.  Although  the  data 
trends  with  frequency  and  amplitude  are  well  captured  by  the 
regression  model,  the  overall  agreement  may  be  compromised 
by  the  fact  that  k  =  0,04  slightly  exceeds  the  expected  range 
of  applicability. 

Four  separate  contributions  were  identified;  i.e.,yni«  =  4.  The 
results  are  summarized  in  Table  3.  As  shown,  the  key  contrib¬ 
utors  are  the  first  and  fourth  parameters,  both  showing  diver¬ 
gence  for  k  ^  0.04.  Thus,  these  results  also  show  that  appli¬ 
cation  of  the  procedure  to  data  taken  at  A:  =  0.04  is  marginal. 
Moreover,  the  regression  analysis  confirms  that  indicial  re¬ 
sponse  time  constants,  as  determined  from  ramp  and  hold  ma¬ 
neuvers,  are  correct  within  approximately  10%.  It  also  suggests 
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that  some  of  the  nonlinear  contributions  are  somewhat  faster, 
contradicting  the  assumption  that  the  decaying  loads  are  simple 
exponentials. 

Conclusions 

Experimental  force  and  moment  measurements  were  made 
for  both  harmonic  and  ramp-and-hold  motions  in  roll.  These 
experiments  were  designed  to  isolate  critical -state  effects.  For 
the  first  time,  evidence  showing  significant  critical-state  tran¬ 
sients  was  presented.  Several  conclusions  regarding  these  crit¬ 
ical-state  transients  are  as  follows. 

1)  Severe  time  lags  ^  30)  existed  for  motions  during 
which  a  critical  state  was  not  crossed  (for  —4  deg  <  (f)  <  4 
deg).  Therefore,  stability  derivatives  are  not  valid  for  reduced 
frequencies,  k  S:  0.033. 

2)  An  independent  assessment  of  the  longest  time-scales  in 
the  1<^|  <  5-deg  range  based  on  relationships  between  stability 
derivatives  and  indicial  response  characteristics  confirmed 

30.  Also,  the  good  agreement  between  independent  ex¬ 
periments  also  supports  the  validity  of  nonlinear  indicial  re¬ 
sponse  theory  as  applied  to  these  data. 

3)  Critical-state  encounters  often  have  a  severe  effect  on  the 
transient  rolling  moment  behavior. 

4)  The  5-deg  critical  state  yielded  the  longest  transient  ef¬ 
fects.  A  relaxation  time,  f  70,  was  found  when  crossing 
this  critical  state  from  either  above  or  below.  Crossing  the 
range  including  the  5-deg  critical  state,  however,  did  riot  yield 
a  noticeable  history  effect  on  the  subsequent  response  for  0 
from  6  to  7  deg. 

5)  The  8-deg  critical  state  yielded  slightly  smaller  relaxation 
times  ^  50  compared  to  the  5-deg  critical  state.  However, 
including  the  flow  history  response  for  a  crossing  of  the  8-deg 
critical  state  yielded  a  noticeable  history  effect  on  the  subse¬ 
quent  response  for  (j>  from  9  to  10  deg. 

6)  The  analysis  of  the  11.3-deg  critical  state  was  not  con¬ 
clusive  because  the  overall  change  in  rolling  moment  coeffi¬ 
cient  from  the  end  of  the  ramp  to  its  steady-state  value  was 
very  small. 

7)  Responses  to  roll  motions  at  the  examined  angle  of  attack 
and  roll  angles  exhibit  three  distinct  time  scales  (potential,  vor¬ 
tical,  and  critical-state  transients).  These  time  scales  must  be 
accounted  for  to  achieve  accurate  aerodynamic  models. 
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